From 6c34f151f6c23ae680727d45755941797972cdb6 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 1 Oct 2026 15:52:00 +0100 Subject: [PATCH 01/33] Implementation of beta fixed point solve, based on Jon M work. --- process/core/caller.py | 112 ++++++++++++++++++++--------------------- 1 file changed, 55 insertions(+), 57 deletions(-) diff --git a/process/core/caller.py b/process/core/caller.py index b2283e1516..87b0316aeb 100644 --- a/process/core/caller.py +++ b/process/core/caller.py @@ -6,6 +6,7 @@ from typing import TYPE_CHECKING import numpy as np +from scipy.optimize import fixed_point from tabulate import tabulate from process.core import constants @@ -18,6 +19,7 @@ from process.data_structure.blanket_variables import BlktModelTypes from process.data_structure.numerics import PROCESSRunMode from process.data_structure.stellarator_variables import StellaratorModel +from process.models.physics.physics import PlasmaBeta from process.models.tfcoil.base import TFConductorModel from process.models.tfcoil.superconducting import SuperconductingTFTurnType @@ -71,67 +73,56 @@ def check_agreement( return np.allclose(previous, current, rtol=1.0e-6, equal_nan=True) return False - def call_models(self, xc: np.ndarray, m: int) -> tuple[float, np.ndarray]: - """Evaluate models until results are idempotent. + def _calculate_beta(self, beta: float, xc: np.ndarray) -> float: + """Return the total plasma beta implied by a model evaluation at ``beta``.""" + self._call_models_once(xc, beta_override=beta) + + beta_thermal = PlasmaBeta.calculate_plasma_beta( + pres_plasma=( + constants.KILOELECTRON_VOLT + * ( + self.data.physics.nd_plasma_electrons_vol_avg + * self.data.physics.temp_plasma_electron_density_weighted_kev + + self.data.physics.nd_plasma_ions_total_vol_avg + * self.data.physics.temp_plasma_ion_density_weighted_kev + ) + ), + b_field=self.data.physics.b_plasma_total, + ) - Ensure objective function and constraints are idempotent before returning. + return ( + self.data.physics.beta_fast_alpha + + self.data.physics.beta_beam + + beta_thermal + ) - Parameters - ---------- - xc : np.ndarray - optimisation parameters - m : int - number of constraints + def _solve_beta_fixed_point(self, xc: np.ndarray) -> None: + """Solve the internally coupled plasma beta fixed-point problem.""" + beta_0 = self.data.physics.beta_total_vol_avg - Returns - ------- - Tuple[float, np.ndarray] - objective function and constraints + try: + beta_solution = fixed_point( + lambda beta: self._calculate_beta(beta, xc), + beta_0, + maxiter=20, + ) + except RuntimeError: + logger.warning("Plasma beta fixed-point solve did not converge; continuing.") + return - Raises - ------ - RuntimeError - if values are non-idempotent after successive - evaluations - """ - objf_prev = None - conf_prev = None - - # Evaluate models up to 10 times; any more implies non-converging values - for _ in range(10): - self._call_models_once(xc) - # Evaluate objective function and constraints - objf = objective_function(self.data.numerics.i_figure_merit, self.data) - conf, _, _, _, _ = constraints.constraint_eqns(m, -1, self.data) - - if objf_prev is None and conf_prev is None: - # First run: run again to check idempotence - logger.debug("New optimisation parameter vector being evaluated") - objf_prev = objf - conf_prev = conf - continue - - # Check for idempotence - if self.check_agreement(objf_prev, objf) and self.check_agreement( - conf_prev, conf - ): - # Idempotent: no longer changing, so return - logger.debug( - "Model evaluations idempotent, returning objective " - "function and constraints" - ) - return objf, conf + # Re-evaluate at the returned fixed point so all model state corresponds + # to the solved beta. + self._calculate_beta(beta_solution, xc) + self.data.physics.beta_total_vol_avg = beta_solution - # Not idempotent: still changing, so evaluate models again - logger.debug("Model evaluations not idempotent: evaluating again") - objf_prev = objf - conf_prev = conf + def call_models(self, xc: np.ndarray, m: int) -> tuple[float, np.ndarray]: + """Evaluate models with plasma beta solved as a fixed-point problem.""" + self._solve_beta_fixed_point(xc) - raise RuntimeError( - "After 10 model evaluations at the current optimisation parameter " - "vector, values for the objective function and constraints haven't " - "converged (don't produce idempotent values)." - ) + objf = objective_function(self.data.numerics.i_figure_merit, self.data) + conf, _, _, _, _ = constraints.constraint_eqns(m, -1, self.data) + + return objf, conf def call_models_and_write_output(self, xc: np.ndarray, ifail: int): """Evaluate models until results are idempotent, then write output files. @@ -164,7 +155,7 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int): # Divert OUT.DAT and MFILE.DAT output to scratch files for # idempotence checking OutputFileManager.open_idempotence_files(self.data.globals.output_prefix) - self._call_models_once(xc) + self._solve_beta_fixed_point(xc) # Write mfile finalise(self.models, self.data, ifail) @@ -247,7 +238,9 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int): non_idempotent_msg=non_idempotent_warning + "\n" + non_idempotent_table, ) - def _call_models_once(self, xc: np.ndarray): + def _call_models_once( + self, xc: np.ndarray, beta_override: float | None = None + ) -> None: """Call the physics and engineering models. This method is the principal caller of all the physics and @@ -268,6 +261,11 @@ def _call_models_once(self, xc: np.ndarray): # Convert variables set_scaled_iteration_variable(xc, nvars, self.data) + # The beta fixed-point solve must override any beta value supplied as an + # iteration variable in xc. + if beta_override is not None: + self.data.physics.beta_total_vol_avg = beta_override + # Perform the various function calls # Stellarator caller if self.data.stellarator.istell != StellaratorModel.DISABLED: From 707589ca1a9c2bd2c1a9edf7b60dd39b6ffe3d39 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Mon, 5 Oct 2026 15:47:17 +0100 Subject: [PATCH 02/33] Revert caller.py to main --- process/core/caller.py | 112 +++++++++++++++++++++-------------------- 1 file changed, 57 insertions(+), 55 deletions(-) diff --git a/process/core/caller.py b/process/core/caller.py index 87b0316aeb..b2283e1516 100644 --- a/process/core/caller.py +++ b/process/core/caller.py @@ -6,7 +6,6 @@ from typing import TYPE_CHECKING import numpy as np -from scipy.optimize import fixed_point from tabulate import tabulate from process.core import constants @@ -19,7 +18,6 @@ from process.data_structure.blanket_variables import BlktModelTypes from process.data_structure.numerics import PROCESSRunMode from process.data_structure.stellarator_variables import StellaratorModel -from process.models.physics.physics import PlasmaBeta from process.models.tfcoil.base import TFConductorModel from process.models.tfcoil.superconducting import SuperconductingTFTurnType @@ -73,56 +71,67 @@ def check_agreement( return np.allclose(previous, current, rtol=1.0e-6, equal_nan=True) return False - def _calculate_beta(self, beta: float, xc: np.ndarray) -> float: - """Return the total plasma beta implied by a model evaluation at ``beta``.""" - self._call_models_once(xc, beta_override=beta) - - beta_thermal = PlasmaBeta.calculate_plasma_beta( - pres_plasma=( - constants.KILOELECTRON_VOLT - * ( - self.data.physics.nd_plasma_electrons_vol_avg - * self.data.physics.temp_plasma_electron_density_weighted_kev - + self.data.physics.nd_plasma_ions_total_vol_avg - * self.data.physics.temp_plasma_ion_density_weighted_kev - ) - ), - b_field=self.data.physics.b_plasma_total, - ) - - return ( - self.data.physics.beta_fast_alpha - + self.data.physics.beta_beam - + beta_thermal - ) + def call_models(self, xc: np.ndarray, m: int) -> tuple[float, np.ndarray]: + """Evaluate models until results are idempotent. - def _solve_beta_fixed_point(self, xc: np.ndarray) -> None: - """Solve the internally coupled plasma beta fixed-point problem.""" - beta_0 = self.data.physics.beta_total_vol_avg + Ensure objective function and constraints are idempotent before returning. - try: - beta_solution = fixed_point( - lambda beta: self._calculate_beta(beta, xc), - beta_0, - maxiter=20, - ) - except RuntimeError: - logger.warning("Plasma beta fixed-point solve did not converge; continuing.") - return + Parameters + ---------- + xc : np.ndarray + optimisation parameters + m : int + number of constraints - # Re-evaluate at the returned fixed point so all model state corresponds - # to the solved beta. - self._calculate_beta(beta_solution, xc) - self.data.physics.beta_total_vol_avg = beta_solution + Returns + ------- + Tuple[float, np.ndarray] + objective function and constraints - def call_models(self, xc: np.ndarray, m: int) -> tuple[float, np.ndarray]: - """Evaluate models with plasma beta solved as a fixed-point problem.""" - self._solve_beta_fixed_point(xc) + Raises + ------ + RuntimeError + if values are non-idempotent after successive + evaluations + """ + objf_prev = None + conf_prev = None + + # Evaluate models up to 10 times; any more implies non-converging values + for _ in range(10): + self._call_models_once(xc) + # Evaluate objective function and constraints + objf = objective_function(self.data.numerics.i_figure_merit, self.data) + conf, _, _, _, _ = constraints.constraint_eqns(m, -1, self.data) + + if objf_prev is None and conf_prev is None: + # First run: run again to check idempotence + logger.debug("New optimisation parameter vector being evaluated") + objf_prev = objf + conf_prev = conf + continue + + # Check for idempotence + if self.check_agreement(objf_prev, objf) and self.check_agreement( + conf_prev, conf + ): + # Idempotent: no longer changing, so return + logger.debug( + "Model evaluations idempotent, returning objective " + "function and constraints" + ) + return objf, conf - objf = objective_function(self.data.numerics.i_figure_merit, self.data) - conf, _, _, _, _ = constraints.constraint_eqns(m, -1, self.data) + # Not idempotent: still changing, so evaluate models again + logger.debug("Model evaluations not idempotent: evaluating again") + objf_prev = objf + conf_prev = conf - return objf, conf + raise RuntimeError( + "After 10 model evaluations at the current optimisation parameter " + "vector, values for the objective function and constraints haven't " + "converged (don't produce idempotent values)." + ) def call_models_and_write_output(self, xc: np.ndarray, ifail: int): """Evaluate models until results are idempotent, then write output files. @@ -155,7 +164,7 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int): # Divert OUT.DAT and MFILE.DAT output to scratch files for # idempotence checking OutputFileManager.open_idempotence_files(self.data.globals.output_prefix) - self._solve_beta_fixed_point(xc) + self._call_models_once(xc) # Write mfile finalise(self.models, self.data, ifail) @@ -238,9 +247,7 @@ def call_models_and_write_output(self, xc: np.ndarray, ifail: int): non_idempotent_msg=non_idempotent_warning + "\n" + non_idempotent_table, ) - def _call_models_once( - self, xc: np.ndarray, beta_override: float | None = None - ) -> None: + def _call_models_once(self, xc: np.ndarray): """Call the physics and engineering models. This method is the principal caller of all the physics and @@ -261,11 +268,6 @@ def _call_models_once( # Convert variables set_scaled_iteration_variable(xc, nvars, self.data) - # The beta fixed-point solve must override any beta value supplied as an - # iteration variable in xc. - if beta_override is not None: - self.data.physics.beta_total_vol_avg = beta_override - # Perform the various function calls # Stellarator caller if self.data.stellarator.istell != StellaratorModel.DISABLED: From 8e6cb331031d92878751e7a8ae1aebc5e6d5ffaf Mon Sep 17 00:00:00 2001 From: ym1906 Date: Mon, 5 Oct 2026 16:03:33 +0100 Subject: [PATCH 03/33] Derive plasma beta explicitly from current plasma state --- process/models/physics/physics.py | 36 ++++++++++++++++++++----------- 1 file changed, 24 insertions(+), 12 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 43975e3f0e..9189f30d20 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -430,12 +430,6 @@ def run(self): # ============================================ - # ----------------------------------------------------- - # Beta Components - # ----------------------------------------------------- - - self.beta.run() - # Stored thermal energies in the plasma ( @@ -694,6 +688,9 @@ def run(self): f_plasma_fuel_deuterium=self.data.physics.f_plasma_fuel_deuterium, ) + # Calculate beta after the fast-particle contributions have been updated. + self.beta.run() + # Calculate ion/electron equilibration power self.data.physics.pden_ion_electron_equilibration_mw = rether( @@ -3506,6 +3503,27 @@ def run(self): If an illegal value of i_beta_norm_max is provided. """ + # Thermal beta is derived directly from the current plasma state. + self.data.physics.beta_thermal_vol_avg = self.calculate_plasma_beta( + pres_plasma=( + constants.KILOELECTRON_VOLT + * ( + self.data.physics.nd_plasma_electrons_vol_avg + * self.data.physics.temp_plasma_electron_density_weighted_kev + + self.data.physics.nd_plasma_ions_total_vol_avg + * self.data.physics.temp_plasma_ion_density_weighted_kev + ) + ), + b_field=self.data.physics.b_plasma_total, + ) + + # Total beta is derived from the thermal and fast-particle components. + self.data.physics.beta_total_vol_avg = ( + self.data.physics.beta_thermal_vol_avg + + self.data.physics.beta_fast_alpha + + self.data.physics.beta_beam + ) + # ----------------------------------------------------- # Normalised Beta Limit # ----------------------------------------------------- @@ -3583,12 +3601,6 @@ def run(self): beta=self.data.physics.beta_total_vol_avg, ) - self.data.physics.beta_thermal_vol_avg = ( - self.data.physics.beta_total_vol_avg - - self.data.physics.beta_fast_alpha - - self.data.physics.beta_beam - ) - self.data.physics.beta_poloidal_eps = ( self.data.physics.beta_poloidal_vol_avg * self.data.physics.eps ) From 8dc8e88a9f08a6870e09083c5ec6cd3cd3ba83c5 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:09:39 +0100 Subject: [PATCH 04/33] Remove constraint 1 as it is not required --- process/core/solver/constraints.py | 38 ------------------------------ 1 file changed, 38 deletions(-) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index c5d786d0f8..f6e3d057cd 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -18,7 +18,6 @@ from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.physics import ( BetaComponentLimits, - PlasmaBeta, PlasmaIgnitionModel, ) from process.models.tfcoil.base import TFConductorModel @@ -216,43 +215,6 @@ def eq(value: float, bound: float, registration: ConstraintRegistration): ) -@ConstraintManager.register_constraint(1, "", "=") -def constraint_equation_1(constraint_registration, data): - """Relationship between beta, temperature (keV) and density - - beta_total_vol_avg: total plasma beta - beta_{ft}: fast alpha beta component - beta_{NBI}: neutral beam beta component - n_e: electron density [/m3] - n_i: total ion density [/m3] - T_e: density weighted average electron temperature [keV] - T_i: density weighted average ion temperature [keV] - B_{tot}: total toroidal + poloidal field [T] - """ - # Density weighted temperature is used here as 〈nT〉 != 〈n〉_V * 〈T〉_V - beta_thermal_total_vol_avg = PlasmaBeta.calculate_plasma_beta( - pres_plasma=( - constants.KILOELECTRON_VOLT - * ( - data.physics.nd_plasma_electrons_vol_avg - * data.physics.temp_plasma_electron_density_weighted_kev - + data.physics.nd_plasma_ions_total_vol_avg - * data.physics.temp_plasma_ion_density_weighted_kev - ) - ), - b_field=data.physics.b_plasma_total, - ) - return eq( - ( - data.physics.beta_fast_alpha - + data.physics.beta_beam - + beta_thermal_total_vol_avg - ), - data.physics.beta_total_vol_avg, - constraint_registration, - ) - - @ConstraintManager.register_constraint(2, "MW/m³", "=") def constraint_equation_2(constraint_registration, data): """ From f01251e9ebb2e774055f1a84f5cd051f750b8410 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:10:14 +0100 Subject: [PATCH 05/33] Remove beta_total_vol_avg as ixc 5 as it is now calculated in physics directly --- process/core/solver/iteration_variables.py | 1 - 1 file changed, 1 deletion(-) diff --git a/process/core/solver/iteration_variables.py b/process/core/solver/iteration_variables.py index fd468c5537..1e4665bb6b 100644 --- a/process/core/solver/iteration_variables.py +++ b/process/core/solver/iteration_variables.py @@ -46,7 +46,6 @@ class IterationVariable: 2: IterationVariable("b_plasma_toroidal_on_axis", "physics", 0.010, 30.00), 3: IterationVariable("rmajor", "physics", 0.1, 50.00), 4: IterationVariable("temp_plasma_electron_vol_avg_kev", "physics", 5.0, 150.0), - 5: IterationVariable("beta_total_vol_avg", "physics", 0.001, 1.0), 6: IterationVariable("nd_plasma_electrons_vol_avg", "physics", 2.0e19, 1.0e21), 7: IterationVariable("f_nd_beam_electron", "physics", 1.0e-6, 1.0), 10: IterationVariable("hfact", "physics", 0.1, 3.0), From ee854790edb3190349ed667718681ec4bb2f95c4 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:10:28 +0100 Subject: [PATCH 06/33] Tidy up removal of constraint 1 --- process/data_structure/numerics.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index 85785576a4..d2d318666e 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -200,7 +200,7 @@ class NumericsData: lablcc: list[str] = field( default_factory=lambda: [ - "⟨β⟩ consistency ", + "NOT USED", "Global power balance consistency ", "Ion power balance ", "Electron power balance ", @@ -295,11 +295,11 @@ class NumericsData: ] ) """Labels describing constraint equations (corresponding itvs)
    - * ( 1) Beta (consistency equation) (itv 5) + * ( 1) NOT USED * ( 2) Global power balance (consistency equation) (itv 10,1,2,3,4,6,11) * ( 3) Ion power balance DEPRECATED (itv 10,1,2,3,4,6,11) * ( 4) Electron power balance DEPRECATED (itv 10,1,2,3,4,6,11) - * ( 5) Density upper limit (itv 9,1,2,3,4,5,6) + * ( 5) Density upper limit (itv 9,1,2,3,4,6) * ( 6) (Epsilon x beta poloidal) upper limit (itv 8,1,2,3,4,6) * ( 7) Beam ion density (NBI) (consistency equation) (itv 7) * ( 8) Neutron wall load upper limit (itv 14,1,2,3,4,6) @@ -409,7 +409,6 @@ class NumericsData: * ( 2) b_plasma_toroidal_on_axis * ( 3) rmajor * ( 4) temp_plasma_electron_vol_avg_kev - * ( 5) beta_total_vol_avg * ( 6) nd_plasma_electrons_vol_avg * ( 7) f_nd_beam_electron * ( 8) NOT USED From 0b1a95bd8e30c65509cbf4012b54d10c2b9d86b5 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:11:02 +0100 Subject: [PATCH 07/33] Remove icc = 1 and ixc = 5 --- tests/regression/input_files/large_tokamak_eval.IN.DAT | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/regression/input_files/large_tokamak_eval.IN.DAT b/tests/regression/input_files/large_tokamak_eval.IN.DAT index 66a44b7c0a..ea66f6dd51 100644 --- a/tests/regression/input_files/large_tokamak_eval.IN.DAT +++ b/tests/regression/input_files/large_tokamak_eval.IN.DAT @@ -3,10 +3,9 @@ i_process_run_mode = -2 *---------------Constraint Equations---------------* * Define number of equality constraints -n_equality_constraints = 2 +n_equality_constraints = 1 * Equalities -icc = 1 * Beta icc = 2 * Global power balance * Inequalities From 6478096ab0386aa1f5ef76eb50cfb6aa827c8251 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:11:07 +0100 Subject: [PATCH 08/33] Remove icc = 1 and ixc = 5 --- tests/regression/input_files/large_tokamak_nof.IN.DAT | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/tests/regression/input_files/large_tokamak_nof.IN.DAT b/tests/regression/input_files/large_tokamak_nof.IN.DAT index 6bdc4128be..8acbca19e2 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 = 2 * Maximum number of solver iterations maxcal = 100 @@ -29,8 +29,6 @@ maxcal = 100 * Beta consistency * *------------------* -icc = 1 -ixc = 5 * beta_total_vol_avg beta_total_vol_avg = 0.03 * Global power balance * From c9362c15c82d73c2f7fb136c8e80561a56b2870a Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:11:12 +0100 Subject: [PATCH 09/33] Remove icc = 1 and ixc = 5 --- tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT b/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT index fc7557a2ca..f103b7046f 100644 --- a/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT +++ b/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT @@ -33,15 +33,12 @@ i_figure_merit = -14 * Switch for figure-of-merit (see lablmm for descriptions * DESCRIPTION: Switch for Figure-of-Merit (1: Minimise major radius) * JUSTIFICATION: Aim to minimise major radius -n_equality_constraints = 4 +n_equality_constraints = 3 * DESCRIPTION: No of equality constraints * JUSTIFICATION: *--------------------------Constraint Equations-------------------------* -icc = 1 -* DESCRIPTION: Beta Consistency -* JUSTIFICATION: Consistency equations should always be on icc = 2 * DESCRIPTION: Global Power Balance Consistency @@ -173,7 +170,6 @@ boundu(4) = 150.0 * DESCRIPTION: Volume averaged electron temperature (keV) * JUSTIFICATION: Iterate temperature to get solution -ixc = 5 beta_total_vol_avg = 3.71965626913446160e-02 * DESCRIPTION: Total Plasma Beta * JUSTIFICATION: Beta is limited by dnbeta From 59315f5e4eabb2aff3aaaf6de4793f5e36d09d70 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:11:16 +0100 Subject: [PATCH 10/33] Remove icc = 1 and ixc = 5 --- tests/regression/input_files/spherical_tokamak_eval.IN.DAT | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT index f79f323f4e..f204a891a0 100644 --- a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT +++ b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT @@ -3,10 +3,9 @@ i_process_run_mode = -2 *---------------Constraint Equations---------------* * Define number of equality constraints -n_equality_constraints = 3 +n_equality_constraints = 2 * Equalities -icc = 1 * Beta icc = 2 * Global power balance icc = 11 * Radial build From 144fe7d005a92f510a73bb0cbebab9f5681b2daa Mon Sep 17 00:00:00 2001 From: ym1906 Date: Tue, 6 Oct 2026 13:11:19 +0100 Subject: [PATCH 11/33] Remove icc = 1 and ixc = 5 --- tests/regression/input_files/st_regression.IN.DAT | 13 +------------ 1 file changed, 1 insertion(+), 12 deletions(-) diff --git a/tests/regression/input_files/st_regression.IN.DAT b/tests/regression/input_files/st_regression.IN.DAT index 5b2f404a4c..a112c0d424 100644 --- a/tests/regression/input_files/st_regression.IN.DAT +++ b/tests/regression/input_files/st_regression.IN.DAT @@ -55,7 +55,7 @@ i_figure_merit = -5 * DESCRIPTION: Switch for Figure-of-Merit * JUSTIFICATION: Selected optimisation -n_equality_constraints = 3 +n_equality_constraints = 2 * DESCRIPTION: Number of equality constraints (the first n_equality_constraints icc = ... statements are considered equality) * JUSTIFICATION: 3 equality (consistency) equations @@ -70,10 +70,6 @@ itart = 1 *_______________________________________________________________________________________* *------------------------------- *Consistency Equations* -------------------------------* *‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾‾* -icc = 1 -* DESCRIPTION: Beta Consistency -* JUSTIFICATION: Consistency equations should always be on -* VARIABLES: Calculated in-situ, no iteration variables icc = 2 * DESCRIPTION: Global Power Balance Consistency @@ -295,13 +291,6 @@ icc = 24 * DESCRIPTION: Plasma beta upper limit * JUSTIFICATION: Turned off, not in use. icc=24. Allow to be calculated in-situ -ixc = 5 -beta_total_vol_avg = 0.2 -boundl(5) = 0.01 -boundu(5) = 1.0 -* DESCRIPTION: Total Plasma Beta, -* JUSTIFICATION: Beta is limited by beta_norm_max - *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ From 122540bf9d0c6b0d10cb0f09fbe7bab2a18b72ce Mon Sep 17 00:00:00 2001 From: ym1906 Date: Wed, 7 Oct 2026 14:19:21 +0100 Subject: [PATCH 12/33] Tidy up beta calculations into plasmabeta.run, and move beta update after fast particle calculation --- process/models/physics/physics.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 9189f30d20..378ef009a4 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -3588,6 +3588,10 @@ def run(self): rminor=self.data.physics.rminor, ) + # ----------------------------------------------------- + # Derived beta components + # ----------------------------------------------------- + self.data.physics.beta_toroidal_vol_avg = ( self.data.physics.beta_total_vol_avg * self.data.physics.b_plasma_total**2 @@ -3613,6 +3617,7 @@ def run(self): ) ** 2 ) + self.data.physics.beta_thermal_toroidal_vol_avg = ( self.data.physics.beta_thermal_vol_avg * ( @@ -3655,6 +3660,7 @@ def run(self): ) ** 2 ) + self.data.physics.beta_norm_poloidal = ( self.data.physics.beta_norm_total * ( From 7cbaa0f7a27821d767a483fe001a2b45bc042a3d Mon Sep 17 00:00:00 2001 From: ym1906 Date: Wed, 7 Oct 2026 14:22:50 +0100 Subject: [PATCH 13/33] Revert "Remove constraint 1 as it is not required" This reverts commit 8dc8e88a9f08a6870e09083c5ec6cd3cd3ba83c5. --- process/core/solver/constraints.py | 38 ++++++++++++++++++++++++++++++ 1 file changed, 38 insertions(+) diff --git a/process/core/solver/constraints.py b/process/core/solver/constraints.py index f6e3d057cd..c5d786d0f8 100644 --- a/process/core/solver/constraints.py +++ b/process/core/solver/constraints.py @@ -18,6 +18,7 @@ from process.models.physics.exhaust import PlasmaExhaust from process.models.physics.physics import ( BetaComponentLimits, + PlasmaBeta, PlasmaIgnitionModel, ) from process.models.tfcoil.base import TFConductorModel @@ -215,6 +216,43 @@ def eq(value: float, bound: float, registration: ConstraintRegistration): ) +@ConstraintManager.register_constraint(1, "", "=") +def constraint_equation_1(constraint_registration, data): + """Relationship between beta, temperature (keV) and density + + beta_total_vol_avg: total plasma beta + beta_{ft}: fast alpha beta component + beta_{NBI}: neutral beam beta component + n_e: electron density [/m3] + n_i: total ion density [/m3] + T_e: density weighted average electron temperature [keV] + T_i: density weighted average ion temperature [keV] + B_{tot}: total toroidal + poloidal field [T] + """ + # Density weighted temperature is used here as 〈nT〉 != 〈n〉_V * 〈T〉_V + beta_thermal_total_vol_avg = PlasmaBeta.calculate_plasma_beta( + pres_plasma=( + constants.KILOELECTRON_VOLT + * ( + data.physics.nd_plasma_electrons_vol_avg + * data.physics.temp_plasma_electron_density_weighted_kev + + data.physics.nd_plasma_ions_total_vol_avg + * data.physics.temp_plasma_ion_density_weighted_kev + ) + ), + b_field=data.physics.b_plasma_total, + ) + return eq( + ( + data.physics.beta_fast_alpha + + data.physics.beta_beam + + beta_thermal_total_vol_avg + ), + data.physics.beta_total_vol_avg, + constraint_registration, + ) + + @ConstraintManager.register_constraint(2, "MW/m³", "=") def constraint_equation_2(constraint_registration, data): """ From 3368eb68603b5f1613cc58ef193e840a23d2af7c Mon Sep 17 00:00:00 2001 From: ym1906 Date: Wed, 7 Oct 2026 14:24:12 +0100 Subject: [PATCH 14/33] Add back in descriptions for icc = 1 --- process/data_structure/numerics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/process/data_structure/numerics.py b/process/data_structure/numerics.py index d2d318666e..4687dddd74 100644 --- a/process/data_structure/numerics.py +++ b/process/data_structure/numerics.py @@ -200,7 +200,7 @@ class NumericsData: lablcc: list[str] = field( default_factory=lambda: [ - "NOT USED", + "⟨β⟩ consistency ", "Global power balance consistency ", "Ion power balance ", "Electron power balance ", @@ -295,7 +295,7 @@ class NumericsData: ] ) """Labels describing constraint equations (corresponding itvs)
      - * ( 1) NOT USED + * ( 1) Beta consistency equation * ( 2) Global power balance (consistency equation) (itv 10,1,2,3,4,6,11) * ( 3) Ion power balance DEPRECATED (itv 10,1,2,3,4,6,11) * ( 4) Electron power balance DEPRECATED (itv 10,1,2,3,4,6,11) From 01d3e3941845108a087de85a7f853f0e57788dfa Mon Sep 17 00:00:00 2001 From: ym1906 Date: Wed, 7 Oct 2026 14:27:32 +0100 Subject: [PATCH 15/33] Remove reference to itv 5 --- process/data_structure/physics_variables.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index dad94478af..ab8c1514b9 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -639,7 +639,7 @@ class PhysicsData: """multiplier for beam-background fusion calculation""" beta_total_vol_avg: float = 0.042 - """Volume averaged total plasma beta (⟨β⟩) (`iteration variable 5`) (calculated if stellarator)""" + """Volume averaged total plasma beta (⟨β⟩) (calculated if stellarator)""" beta_fast_alpha: float = 0.0 """Fast alpha beta component (β_alpha)""" From b90996c2a2c257ac3dfae20357e6d22b754f8a23 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Wed, 7 Oct 2026 14:30:28 +0100 Subject: [PATCH 16/33] Remove references to plasma beta consistency equations --- documentation/source/solver/solver-guide.md | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/documentation/source/solver/solver-guide.md b/documentation/source/solver/solver-guide.md index be891cce1a..82e5316bf5 100644 --- a/documentation/source/solver/solver-guide.md +++ b/documentation/source/solver/solver-guide.md @@ -18,7 +18,6 @@ A PROCESS input file will, for example, define which constraint equations are be n_equality_constraints = 3 * Equalities -icc = 1 * Beta icc = 2 * Global power balance icc = 11 * Radial build @@ -144,7 +143,6 @@ i_process_run_mode = -2 * evaluation mode n_equality_constraints = 2 * Equalities -icc = 1 * Beta icc = 2 * Global power balance * Inequalities @@ -160,4 +158,4 @@ ixc = 4 * temp_plasma_electron_vol_avg_kev ixc = 6 * nd_plasma_electrons_vol_avg ... ``` -The beta and plasma power balance equality constraints are recommended to ensure model consistency, and the `temp_plasma_electron_vol_avg_kev` and `nd_plasma_electrons_vol_avg` iteration variables used to solve them. +The plasma power balance equality constraint is recommended to ensure model consistency, and the `temp_plasma_electron_vol_avg_kev` and `nd_plasma_electrons_vol_avg` may be used as iteration variables used to solve them. From a6bc5b05264831f398b08167ced8c693bbd48f47 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Wed, 7 Oct 2026 14:41:32 +0100 Subject: [PATCH 17/33] Update plasma_beta docs to reflect new status of icc = 1 --- .../source/physics-models/plasma_beta/plasma_beta.md | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/documentation/source/physics-models/plasma_beta/plasma_beta.md b/documentation/source/physics-models/plasma_beta/plasma_beta.md index 20e81f5a14..6eafc171ad 100644 --- a/documentation/source/physics-models/plasma_beta/plasma_beta.md +++ b/documentation/source/physics-models/plasma_beta/plasma_beta.md @@ -390,7 +390,7 @@ $$ ### Beta consistency -This constraint can be activated by stating `icc = 1` in the input file. +beta_total_vol_avg is calculated directly from the current plasma state and is no longer an iteration variable. Constraint 1 may be activated as an optional global consistency check or safeguard against potential non-idempotent behaviour by stating icc = 1 in the input file. Ensures the relationship between $\beta$, density, temperature and total magnetic field is withheld by checking the fixed input or iteration variable $\texttt{beta}$ is consistent in value with the rest of the physics parameters @@ -404,8 +404,6 @@ $$ \langle n_{\text{e}}T_{\text{e}} \rangle_{\text{V}} = \overbrace{\langle n_{\text{e}} \rangle_{\text{V}}}^{\texttt{nd_plasma_electrons_vol_avg}} \times \overbrace{\langle T_{\text{e}} \rangle_{\text{n}}}^{\texttt{temp_plasma_electron_density_weighted_kev}} $$ -**It is highly recommended to always have this constraint on as it is a global consistency checker** - ---------------- ### Poloidal beta and inverse aspect upper limit From d449469fb9b1559e630807ee69e872b446f87f97 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Wed, 7 Oct 2026 14:41:57 +0100 Subject: [PATCH 18/33] Tidy up references to icc = 1 and ixc = 5 --- documentation/source/physics-models/profiles/plasma_profiles.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/source/physics-models/profiles/plasma_profiles.md b/documentation/source/physics-models/profiles/plasma_profiles.md index 362da0f76e..4697d4f4df 100644 --- a/documentation/source/physics-models/profiles/plasma_profiles.md +++ b/documentation/source/physics-models/profiles/plasma_profiles.md @@ -405,7 +405,7 @@ $$ With the coefficients used to turn the temperature from $\text{keV}$ back to Joules. -We calculate the volume averaged thermal plasma pressure using the volume averaged density and the density weighted volume averaged temperature for the electrons and ions. The density weighted temperature is used as $\langle nT \rangle \neq \langle n \rangle \langle T \rangle$. This can be seen also in the [equality constraint for thermal $\beta$](../plasma_beta/plasma_beta.md#beta-consistency) +We calculate the volume averaged thermal plasma pressure using the volume averaged density and the density weighted volume averaged temperature for the electrons and ions. The density weighted temperature is used as $\langle nT \rangle \neq \langle n \rangle \langle T \rangle$. $$ \langle p_{\text{thermal}} \rangle = \overbrace{\langle n_{\text{e}} \rangle_{\text{V}}}^{\texttt{nd_plasma_electrons_vol_avg}} \times \overbrace{\langle T_{\text{e}} \rangle_{\text{n}}}^{\texttt{temp_plasma_electron_density_weighted_kev}}\\ From ed7e2afb1db0bebc23399cf12826c4c7a3fe9119 Mon Sep 17 00:00:00 2001 From: Graeme Turkington <107113942+grmtrkngtn@users.noreply.github.com> Date: Thu, 8 Oct 2026 11:02:11 +0100 Subject: [PATCH 19/33] Update process/data_structure/physics_variables.py Co-authored-by: Timothy <75321887+timothy-nunn@users.noreply.github.com> --- process/data_structure/physics_variables.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index ab8c1514b9..bbc0b8573d 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -639,7 +639,7 @@ class PhysicsData: """multiplier for beam-background fusion calculation""" beta_total_vol_avg: float = 0.042 - """Volume averaged total plasma beta (⟨β⟩) (calculated if stellarator)""" + """Volume averaged total plasma beta (⟨β⟩)""" beta_fast_alpha: float = 0.0 """Fast alpha beta component (β_alpha)""" From c82ac5cf17122419d7fe8f8972983a33b7070b5d Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 11:00:30 +0100 Subject: [PATCH 20/33] Tidied up the end of solver guide to tell user to select one ixc --- documentation/source/solver/solver-guide.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/source/solver/solver-guide.md b/documentation/source/solver/solver-guide.md index 82e5316bf5..a4d0ecfebe 100644 --- a/documentation/source/solver/solver-guide.md +++ b/documentation/source/solver/solver-guide.md @@ -158,4 +158,4 @@ ixc = 4 * temp_plasma_electron_vol_avg_kev ixc = 6 * nd_plasma_electrons_vol_avg ... ``` -The plasma power balance equality constraint is recommended to ensure model consistency, and the `temp_plasma_electron_vol_avg_kev` and `nd_plasma_electrons_vol_avg` may be used as iteration variables used to solve them. +The plasma power balance equality constraint is recommended to ensure model consistency, the user may select one of `temp_plasma_electron_vol_avg_kev` and `nd_plasma_electrons_vol_avg` to solve it. From cfea93ed5e33528869640673c6cd76140be8bc0f Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 11:03:24 +0100 Subject: [PATCH 21/33] Removed beta_total_vol_avg as input --- process/core/input.py | 1 - 1 file changed, 1 deletion(-) diff --git a/process/core/input.py b/process/core/input.py index 31b6c6f6cd..73f67ac078 100644 --- a/process/core/input.py +++ b/process/core/input.py @@ -164,7 +164,6 @@ def bounds(self) -> tuple[NumberType | None, NumberType | None]: "alphat": InputVariable("physics", float, range=(0.0, 10.0)), "aspect": InputVariable("physics", float, range=(1.001, 40.0)), "beamfus0": InputVariable("physics", float, range=(0.01, 10.0)), - "beta_total_vol_avg": InputVariable("physics", float, range=(0.0, 1.0)), "beta_vol_avg_max": InputVariable("physics", float, range=(0.0, 1.0)), "beta_vol_avg_min": InputVariable("physics", float, range=(0.0, 1.0)), "betbm0": InputVariable("physics", float, range=(0.0, 10.0)), From caa8643229f66f9a34d37cf41b3fa3a3db7b8818 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 11:14:21 +0100 Subject: [PATCH 22/33] Removed beta_total_vol_avg as input, removed density as ixc as only 1 equality contraint --- tests/regression/input_files/large_tokamak_eval.IN.DAT | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/regression/input_files/large_tokamak_eval.IN.DAT b/tests/regression/input_files/large_tokamak_eval.IN.DAT index ea66f6dd51..e429ad9720 100644 --- a/tests/regression/input_files/large_tokamak_eval.IN.DAT +++ b/tests/regression/input_files/large_tokamak_eval.IN.DAT @@ -37,7 +37,7 @@ icc = 32 * TF coil conduit stress upper limit ixc = 4 * temp_plasma_electron_vol_avg_kev boundu(4) = 100.0 -ixc = 6 * nd_plasma_electrons_vol_avg + *---------------Cs Fatigue Variables---------------* @@ -273,7 +273,6 @@ alphan = 1.00 * density profile index alphat = 1.45 * temperature profile index i_alphaj = 1 * Current profile index selection (0 = user input) aspect = 3.0 * aspect ratio (`iteration variable 1`) -beta_total_vol_avg = 0.03230408815355488 * total plasma beta (`iteration variable 5`) (calculated if stellarator) b_plasma_toroidal_on_axis = 5.318322174644904 * toroidal field on axis (T) (`iteration variable 2`) nd_plasma_electrons_vol_avg = 7.796223900029837e+19 * electron density (/m3) (`iteration variable 6`) beta_norm_max = 3.0 * Troyon-like coefficient for beta scaling calculated From 779034db9eb5a7d923709792d245aa20c1e44111 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 11:14:51 +0100 Subject: [PATCH 23/33] Removed beta_total_vol_avg as input file --- tests/regression/input_files/large_tokamak_nof.IN.DAT | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/regression/input_files/large_tokamak_nof.IN.DAT b/tests/regression/input_files/large_tokamak_nof.IN.DAT index 8acbca19e2..ac584697b6 100644 --- a/tests/regression/input_files/large_tokamak_nof.IN.DAT +++ b/tests/regression/input_files/large_tokamak_nof.IN.DAT @@ -27,10 +27,6 @@ maxcal = 100 * Constraint Equations - Consistency Equations * ************************************************ -* Beta consistency * -*------------------* -beta_total_vol_avg = 0.03 - * Global power balance * *----------------------* icc = 2 From 6a21140f736f2bf57051d85f3ee9dabe117d142b Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 11:15:16 +0100 Subject: [PATCH 24/33] Removed beta_total_vol_avg as input --- tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT | 4 ---- 1 file changed, 4 deletions(-) diff --git a/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT b/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT index f103b7046f..f58ad98972 100644 --- a/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT +++ b/tests/regression/input_files/low_aspect_ratio_DEMO.IN.DAT @@ -170,10 +170,6 @@ boundu(4) = 150.0 * DESCRIPTION: Volume averaged electron temperature (keV) * JUSTIFICATION: Iterate temperature to get solution -beta_total_vol_avg = 3.71965626913446160e-02 -* DESCRIPTION: Total Plasma Beta -* JUSTIFICATION: Beta is limited by dnbeta - ixc = 6 nd_plasma_electrons_vol_avg = 6.88360041658364314e+19 * DESCRIPTION: Electron density (/m3) From da869379df5d1a4489d65ff2067ba82613786b82 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 11:15:48 +0100 Subject: [PATCH 25/33] removed electron density as ixc to correctly form system of equations --- tests/regression/input_files/spherical_tokamak_eval.IN.DAT | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT index f204a891a0..a7a8282f58 100644 --- a/tests/regression/input_files/spherical_tokamak_eval.IN.DAT +++ b/tests/regression/input_files/spherical_tokamak_eval.IN.DAT @@ -31,9 +31,7 @@ icc = 16 * Net electric lower limit ixc = 4 * temp_plasma_electron_vol_avg_kev boundl(4) = 5.0 boundu(4) = 25.0 -ixc = 6 * nd_plasma_electrons_vol_avg -boundl(6) = 0.5d+20 -boundu(6) = 5.0d+20 + ixc = 29 * bore boundl(29) = 0.1 boundu(29) = 0.8 @@ -257,7 +255,6 @@ i_alphaj = 0 * Current profile index selection (0 = user input) alphan = 0.9 * density profile index alphat = 1.4 * temperature profile index aspect = 1.8 * aspect ratio (`iteration variable 1`) -beta_total_vol_avg = 0.13134204235647895 * total plasma beta (`iteration variable 5`) (calculated if stellarator) b_plasma_toroidal_on_axis = 3.0 * toroidal field on axis (T) (`iteration variable 2`) nd_plasma_electrons_vol_avg = 9.69888313737236e+19 * electron density (/m3) (`iteration variable 6`) beta_norm_max = 5.0 * Troyon-like coefficient for beta scaling From b8b906a84a4c28fe5bbf732b07c1b8f1fe734359 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 13:54:44 +0100 Subject: [PATCH 26/33] copy latest large tokamak eval from regression test over --- examples/data/large_tokamak_eval_IN.DAT | 18 ++++++++++-------- 1 file changed, 10 insertions(+), 8 deletions(-) diff --git a/examples/data/large_tokamak_eval_IN.DAT b/examples/data/large_tokamak_eval_IN.DAT index 96fe08dcc4..e429ad9720 100644 --- a/examples/data/large_tokamak_eval_IN.DAT +++ b/examples/data/large_tokamak_eval_IN.DAT @@ -3,10 +3,9 @@ i_process_run_mode = -2 *---------------Constraint Equations---------------* * Define number of equality constraints -n_equality_constraints = 2 +n_equality_constraints = 1 * Equalities -icc = 1 * Beta icc = 2 * Global power balance * Inequalities @@ -38,7 +37,7 @@ icc = 32 * TF coil conduit stress upper limit ixc = 4 * temp_plasma_electron_vol_avg_kev boundu(4) = 100.0 -ixc = 6 * nd_plasma_electrons_vol_avg + *---------------Cs Fatigue Variables---------------* @@ -101,7 +100,7 @@ p_fusion_total_max_mw = 3000 * maximum fusion power (MW) (`constraint equation p_div_bt_q_aspect_rmajor_max_mw = 10.0 * maximum ratio of Psep*Bt/qAR (MWT/m) (`constraint equation 68`) t_burn_min = 7200.0 * minimum burn time (s) (KE - no longer itv;; see issue #706) pflux_fw_neutron_max_mw = 2.0 * allowable neutron wall-load (MW/m2) (`constraint equation 8`) -f_t_alpha_energy_confinement_min = 5.0 * Lower limit on f_t_alpha_energy_confinementnt the ratio of alpha particle to energy confinement +f_t_alpha_energy_confinement_min = 5.0 * Lower limit on f_t_alpha_energy_confinement the ratio of alpha particle to energy confinement *-------------------Constraints--------------------* @@ -272,8 +271,8 @@ zref(10) = 1.0 alphan = 1.00 * density profile index alphat = 1.45 * temperature profile index +i_alphaj = 1 * Current profile index selection (0 = user input) aspect = 3.0 * aspect ratio (`iteration variable 1`) -beta_total_vol_avg = 0.03230408815355488 * total plasma beta (`iteration variable 5`) (calculated if stellarator) b_plasma_toroidal_on_axis = 5.318322174644904 * toroidal field on axis (T) (`iteration variable 2`) nd_plasma_electrons_vol_avg = 7.796223900029837e+19 * electron density (/m3) (`iteration variable 6`) beta_norm_max = 3.0 * Troyon-like coefficient for beta scaling calculated @@ -294,9 +293,9 @@ plasma_res_factor = 0.7 * plasma resistivity pre-factor radius_plasma_pedestal_density_norm = 0.94 * r/a of density pedestal (`i_plasma_pedestal==1`) radius_plasma_pedestal_temp_norm = 0.94 * r/a of temperature pedestal (`i_plasma_pedestal==1`) tbeta = 2.0 * temperature profile index beta (`i_plasma_pedestal==1) -temp_plasma_separatrix_electron_kev = 5.5 * electron temperature of pedestal (keV) (`i_plasma_pedestal==1') +temp_plasma_pedestal_electron_kev = 5.5 * electron temperature of pedestal (keV) (`i_plasma_pedestal==1') temp_plasma_separatrix_electron_kev = 0.1 * electron temperature at separatrix (keV) (`i_plasma_pedestal==1`) calculated if reinke -i_confinement_time = 34 * switch for energy confinement time scaling law (see description in `tauscl`) +i_confinement_time = 34 * switch for energy confinement time scaling law (see description in `labels_confinement_scalings`) i_plasma_geometry = 0 * switch for plasma cross-sectional shape calculation; kappa = 1.85 * plasma separatrix elongation (calculated if `i_plasma_geometry = 1-5; 7 or 9-10`) q95 = 3.7339078193128556 * Safety factor 'near' plasma edge (`iteration variable 18`) equal to q95 @@ -307,6 +306,7 @@ i_single_null = 1 * switch for single null / double null plasma; f_sync_reflect = 0.6 * synchrotron wall reflectivity factor temp_plasma_electron_vol_avg_kev = 12.221383528378944 * volume averaged electron temperature (keV) (`iteration variable 4`) 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) *----------------------Power-----------------------* @@ -375,13 +375,15 @@ t_tf_superconductor_quench = 17.97282589344206 * fast discharge time for TF co n_tf_coils = 16 * Number of TF coils (default = 50 for stellarators); Number of TF coils outer legs for ST tftmp = 4.75 * peak helium coolant temperature in TF coils and PF coils (K) dr_tf_nose_case = 0.2816873221155309 * inboard TF coil case outer (non-plasma side) thickness (m) (`iteration variable 57`) -dr_tf_wp_with_insulation = 0.5153787768966674 * radial thickness of winding pack (m) (`iteration variable 140`) (issue #514) +dr_tf_wp_with_insulation = 0.5153787768966674 * radial thickness of winding pack (m) (issue #514) dx_tf_turn_steel = 0.008012110032981922 * TF coil conduit case thickness (m) (`iteration variable 58`) dx_tf_wp_insulation = 0.008 * Thickness of the ground insulation layer surrounding (m) temp_cs_superconductor_margin_min = 1.5 * minimum allowable temperature margin ; CS (K) tmargmin = 1.5 * minimum allowable temperature margin ; TFC AND CS (K) v_tf_coil_dump_quench_max_kv = 10.0 * max voltage across TF coil during quench (kV) (`iteration variable 52`) f_a_tf_turn_cable_space_extra_void = 0.3 * coolant fraction of TFC 'cable' (`i_tf_sup=1`); or of TFC leg (`i_tf_sup=0`) +flu_tf_neutron_fast_max = 1e22 * max fluence in the TF coil + *-----------------Times Variables------------------* From be3e51807a03958fb2590c833fe90419f4b8cc0e Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 13:55:31 +0100 Subject: [PATCH 27/33] remove icc = 1 and ixc 5 and beta_total_vol_avg --- examples/data/large_tokamak_varyrun_IN.DAT | 7 +------ 1 file changed, 1 insertion(+), 6 deletions(-) diff --git a/examples/data/large_tokamak_varyrun_IN.DAT b/examples/data/large_tokamak_varyrun_IN.DAT index 1b76e61add..567e4a2e79 100644 --- a/examples/data/large_tokamak_varyrun_IN.DAT +++ b/examples/data/large_tokamak_varyrun_IN.DAT @@ -17,16 +17,11 @@ i_figure_merit = 1 * Error tolerance for VMCON epsvmc = 1e-7 -n_equality_constraints = 3 +n_equality_constraints = 2 * Constraint Equations - Consistency Equations * ************************************************ -* Beta consistency * -*------------------* -icc = 1 -ixc = 5 * beta_total_vol_avg -beta_total_vol_avg = 0.03 * Global power balance * *----------------------* From 61bf24a905917fe5c0898f6d275ac27d66edba6b Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 13:55:54 +0100 Subject: [PATCH 28/33] remove icc = 1 and ixc 5 and beta_total_vol_avg --- .../data/large_tokamak_varied_min_net_electric_IN.DAT | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/examples/data/large_tokamak_varied_min_net_electric_IN.DAT b/examples/data/large_tokamak_varied_min_net_electric_IN.DAT index f60f49cdb1..652aea9d3e 100644 --- a/examples/data/large_tokamak_varied_min_net_electric_IN.DAT +++ b/examples/data/large_tokamak_varied_min_net_electric_IN.DAT @@ -17,17 +17,11 @@ i_figure_merit = 1 * Error tolerance for VMCON epsvmc = 1e-7 -n_equality_constraints = 3 +n_equality_constraints = 1 * Constraint Equations - Consistency Equations * ************************************************ -* Beta consistency * -*------------------* -icc = 1 -ixc = 5 * beta_total_vol_avg -beta_total_vol_avg = 0.03 - * Global power balance * *----------------------* icc = 2 From 79917aa8b184c553585fcc8fa6cfbabccb87efd1 Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 13:57:51 +0100 Subject: [PATCH 29/33] remove icc = 1 and ixc 5 and beta_total_vol_avg --- tests/integration/data/large_tokamak_eval.IN.DAT | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/tests/integration/data/large_tokamak_eval.IN.DAT b/tests/integration/data/large_tokamak_eval.IN.DAT index 852633fa65..25581eb529 100644 --- a/tests/integration/data/large_tokamak_eval.IN.DAT +++ b/tests/integration/data/large_tokamak_eval.IN.DAT @@ -3,10 +3,9 @@ i_process_run_mode = -2 *---------------Constraint Equations---------------* * Define number of equality constraints -n_equality_constraints = 2 +n_equality_constraints = 1 * Equalities -icc = 1 * Beta icc = 2 * Global power balance * Inequalities From 55b1f4d19df032b0e37a432a79a20503067a88c9 Mon Sep 17 00:00:00 2001 From: Graeme Turkington <107113942+grmtrkngtn@users.noreply.github.com> Date: Thu, 8 Oct 2026 15:40:50 +0100 Subject: [PATCH 30/33] Update documentation/source/physics-models/plasma_beta/plasma_beta.md Co-authored-by: Christopher Ashe <91618944+chris-ashe@users.noreply.github.com> --- documentation/source/physics-models/plasma_beta/plasma_beta.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/source/physics-models/plasma_beta/plasma_beta.md b/documentation/source/physics-models/plasma_beta/plasma_beta.md index 6eafc171ad..634adc1850 100644 --- a/documentation/source/physics-models/plasma_beta/plasma_beta.md +++ b/documentation/source/physics-models/plasma_beta/plasma_beta.md @@ -390,7 +390,7 @@ $$ ### Beta consistency -beta_total_vol_avg is calculated directly from the current plasma state and is no longer an iteration variable. Constraint 1 may be activated as an optional global consistency check or safeguard against potential non-idempotent behaviour by stating icc = 1 in the input file. +`beta_total_vol_avg` is calculated directly from the current plasma state and is no longer an iteration variable. Constraint 1 may be activated as an optional global consistency check or safeguard against potential non-idempotent behaviour by stating `icc = 1` in the input file. Ensures the relationship between $\beta$, density, temperature and total magnetic field is withheld by checking the fixed input or iteration variable $\texttt{beta}$ is consistent in value with the rest of the physics parameters From 4508b31463566be802acf593ae120b879dce33dc Mon Sep 17 00:00:00 2001 From: ym1906 Date: Thu, 8 Oct 2026 16:26:11 +0100 Subject: [PATCH 31/33] Moved fast alpha beta calculation into beta class --- process/models/physics/physics.py | 33 +++++++++++++++++-------------- 1 file changed, 18 insertions(+), 15 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 378ef009a4..1a442b8ba4 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -674,21 +674,7 @@ def run(self): self.data.physics.f_p_alpha_plasma_deposited, ) - self.data.physics.beta_fast_alpha = self.beta.fast_alpha_beta( - b_plasma_poloidal_average=self.data.physics.b_plasma_surface_poloidal_average, - b_plasma_toroidal_on_axis=self.data.physics.b_plasma_toroidal_on_axis, - nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg, - nd_plasma_fuel_ions_vol_avg=self.data.physics.nd_plasma_fuel_ions_vol_avg, - nd_plasma_ions_total_vol_avg=self.data.physics.nd_plasma_ions_total_vol_avg, - temp_plasma_electron_density_weighted_kev=self.data.physics.temp_plasma_electron_density_weighted_kev, - temp_plasma_ion_density_weighted_kev=self.data.physics.temp_plasma_ion_density_weighted_kev, - pden_alpha_total_vol_avg_mw=self.data.physics.pden_alpha_total_vol_avg_mw, - pden_plasma_alpha_vol_avg_mw=self.data.physics.pden_plasma_alpha_vol_avg_mw, - i_beta_fast_alpha=self.data.physics.i_beta_fast_alpha, - f_plasma_fuel_deuterium=self.data.physics.f_plasma_fuel_deuterium, - ) - - # Calculate beta after the fast-particle contributions have been updated. + # Calculate beta quantities. self.beta.run() # Calculate ion/electron equilibration power @@ -3517,6 +3503,23 @@ def run(self): b_field=self.data.physics.b_plasma_total, ) + # Fast-alpha beta is calculated here so all tokamak beta calculations are + # owned by PlasmaBeta.run(). Keep fast_alpha_beta() unchanged because it is + # also called directly by the stellarator physics model. + self.data.physics.beta_fast_alpha = self.fast_alpha_beta( + b_plasma_poloidal_average=self.data.physics.b_plasma_surface_poloidal_average, + b_plasma_toroidal_on_axis=self.data.physics.b_plasma_toroidal_on_axis, + nd_plasma_electrons_vol_avg=self.data.physics.nd_plasma_electrons_vol_avg, + nd_plasma_fuel_ions_vol_avg=self.data.physics.nd_plasma_fuel_ions_vol_avg, + nd_plasma_ions_total_vol_avg=self.data.physics.nd_plasma_ions_total_vol_avg, + temp_plasma_electron_density_weighted_kev=self.data.physics.temp_plasma_electron_density_weighted_kev, + temp_plasma_ion_density_weighted_kev=self.data.physics.temp_plasma_ion_density_weighted_kev, + pden_alpha_total_vol_avg_mw=self.data.physics.pden_alpha_total_vol_avg_mw, + pden_plasma_alpha_vol_avg_mw=self.data.physics.pden_plasma_alpha_vol_avg_mw, + i_beta_fast_alpha=self.data.physics.i_beta_fast_alpha, + f_plasma_fuel_deuterium=self.data.physics.f_plasma_fuel_deuterium, + ) + # Total beta is derived from the thermal and fast-particle components. self.data.physics.beta_total_vol_avg = ( self.data.physics.beta_thermal_vol_avg From af9fdcec25d7802dd585b51af78ccde2b74e46bd Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 14 Oct 2026 11:25:32 +0100 Subject: [PATCH 32/33] Remove ixc=5 and beta_total_vol_avg from test input files --- examples/data/large_tokamak_IN.DAT | 8 +------- examples/data/scan_example_file_IN.DAT | 8 +------- tests/integration/data/large_tokamak_IN.DAT | 8 +------- tests/integration/data/large_tokamak_eval.IN.DAT | 1 - tests/unit/data/large_tokamak_IN.DAT | 8 +------- 5 files changed, 4 insertions(+), 29 deletions(-) diff --git a/examples/data/large_tokamak_IN.DAT b/examples/data/large_tokamak_IN.DAT index 00f98e25e7..cc19dd2588 100644 --- a/examples/data/large_tokamak_IN.DAT +++ b/examples/data/large_tokamak_IN.DAT @@ -17,17 +17,11 @@ i_figure_merit = 1 * Error tolerance for VMCON epsvmc = 1e-7 -n_equality_constraints = 3 +n_equality_constraints = 2 * Constraint Equations - Consistency Equations * ************************************************ -* Beta consistency * -*------------------* -icc = 1 -ixc = 5 * beta_total_vol_avg -beta_total_vol_avg = 0.03 - * Global power balance * *----------------------* icc = 2 diff --git a/examples/data/scan_example_file_IN.DAT b/examples/data/scan_example_file_IN.DAT index c93e031720..323c0321ca 100644 --- a/examples/data/scan_example_file_IN.DAT +++ b/examples/data/scan_example_file_IN.DAT @@ -20,7 +20,7 @@ epsvmc = 1e-7 * Inequality constraint tolerance force_vmcon_inequality_tolerance = 1e-6 -n_equality_constraints = 3 +n_equality_constraints = 2 * Scan Variables * ****************** @@ -32,12 +32,6 @@ sweep = 10.5, 10.4, 10.3, 10.2, 10.1, 10.0 * Constraint Equations - Consistency Equations * ************************************************ -* Beta consistency * -*------------------* -icc = 1 -ixc = 5 * beta_total_vol_avg -beta_total_vol_avg = 0.03 - * Global power balance * *----------------------* icc = 2 diff --git a/tests/integration/data/large_tokamak_IN.DAT b/tests/integration/data/large_tokamak_IN.DAT index fe77b15abf..16230592e5 100644 --- a/tests/integration/data/large_tokamak_IN.DAT +++ b/tests/integration/data/large_tokamak_IN.DAT @@ -17,17 +17,11 @@ i_figure_merit = 1 * Error tolerance for VMCON epsvmc = 1e-7 -n_equality_constraints = 3 +n_equality_constraints = 2 * Constraint Equations - Consistency Equations * ************************************************ -* Beta consistency * -*------------------* -icc = 1 -ixc = 5 * beta_total_vol_avg -beta_total_vol_avg = 0.03 - * Global power balance * *----------------------* icc = 2 diff --git a/tests/integration/data/large_tokamak_eval.IN.DAT b/tests/integration/data/large_tokamak_eval.IN.DAT index 25581eb529..46621e030d 100644 --- a/tests/integration/data/large_tokamak_eval.IN.DAT +++ b/tests/integration/data/large_tokamak_eval.IN.DAT @@ -272,7 +272,6 @@ zref(10) = 1.0 alphan = 1.00 * density profile index alphat = 1.45 * temperature profile index aspect = 3.0 * aspect ratio (`iteration variable 1`) -beta_total_vol_avg = 0.03230408815355488 * total plasma beta (`iteration variable 5`) (calculated if stellarator) b_plasma_toroidal_on_axis = 5.318322174644904 * toroidal field on axis (T) (`iteration variable 2`) nd_plasma_electrons_vol_avg = 7.796223900029837e+19 * electron density (/m3) (`iteration variable 6`) beta_norm_max = 3.0 * Troyon-like coefficient for beta scaling calculated diff --git a/tests/unit/data/large_tokamak_IN.DAT b/tests/unit/data/large_tokamak_IN.DAT index fe77b15abf..16230592e5 100644 --- a/tests/unit/data/large_tokamak_IN.DAT +++ b/tests/unit/data/large_tokamak_IN.DAT @@ -17,17 +17,11 @@ i_figure_merit = 1 * Error tolerance for VMCON epsvmc = 1e-7 -n_equality_constraints = 3 +n_equality_constraints = 2 * Constraint Equations - Consistency Equations * ************************************************ -* Beta consistency * -*------------------* -icc = 1 -ixc = 5 * beta_total_vol_avg -beta_total_vol_avg = 0.03 - * Global power balance * *----------------------* icc = 2 From 7dcf9c3868c29d6b87d10a2bb9294cf7985c631a Mon Sep 17 00:00:00 2001 From: Timothy Nunn Date: Wed, 14 Oct 2026 11:49:56 +0100 Subject: [PATCH 33/33] Temporarily remove density iteration from integration eval file --- tests/integration/data/large_tokamak_eval.IN.DAT | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/integration/data/large_tokamak_eval.IN.DAT b/tests/integration/data/large_tokamak_eval.IN.DAT index 46621e030d..b739105f1f 100644 --- a/tests/integration/data/large_tokamak_eval.IN.DAT +++ b/tests/integration/data/large_tokamak_eval.IN.DAT @@ -37,7 +37,7 @@ icc = 32 * TF coil conduit stress upper limit ixc = 4 * temp_plasma_electron_vol_avg_kev boundu(4) = 100.0 -ixc = 6 * nd_plasma_electrons_vol_avg +*ixc = 6 * nd_plasma_electrons_vol_avg *---------------Cs Fatigue Variables---------------*