Derive plasma beta explicitly from the current plasma state

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

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}}\\

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, the user may select one of `temp_plasma_electron_vol_avg_kev` and `nd_plasma_electrons_vol_avg` to solve it.

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

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

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

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

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

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)),

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),

process/data_structure/numerics.py

@@ -295,11 +295,11 @@ class NumericsData:
         ]
     )
     """Labels describing constraint equations (corresponding itvs)<UL>
-    * ( 1) Beta (consistency equation) (itv 5)
+    * ( 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)
-    * ( 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

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 (⟨β⟩)"""
 
     beta_fast_alpha: float = 0.0
     """Fast alpha beta component (β_alpha)"""

process/models/physics/physics.py

@@ -430,12 +430,6 @@ def run(self):
 
         # ============================================
 
-        # -----------------------------------------------------
-        # Beta Components
-        # -----------------------------------------------------
-
-        self.beta.run()
-
         # Stored thermal energies in the plasma
 
         (
@@ -680,19 +674,8 @@ 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 quantities.
+        self.beta.run()
 
         # Calculate ion/electron equilibration power
 
@@ -3506,6 +3489,44 @@ 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,
+        )
+
+        # 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
+            + self.data.physics.beta_fast_alpha
+            + self.data.physics.beta_beam
+        )
+
         # -----------------------------------------------------
         # Normalised Beta Limit
         # -----------------------------------------------------
@@ -3570,6 +3591,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
@@ -3583,12 +3608,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
         )
@@ -3601,6 +3620,7 @@ def run(self):
             )
             ** 2
         )
+
         self.data.physics.beta_thermal_toroidal_vol_avg = (
             self.data.physics.beta_thermal_vol_avg
             * (
@@ -3643,6 +3663,7 @@ def run(self):
             )
             ** 2
         )
+
         self.data.physics.beta_norm_poloidal = (
             self.data.physics.beta_norm_total
             * (

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