From 7f695be13ec5b0709ebf5ba66982d5308acc7f8b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 13:52:59 +0000 Subject: [PATCH 01/23] Add PlasmaBeta class and integrate into Physics initialization --- process/main.py | 7 +++++-- process/physics.py | 11 ++++++++++- 2 files changed, 15 insertions(+), 3 deletions(-) diff --git a/process/main.py b/process/main.py index 029aee551..4d86437ee 100644 --- a/process/main.py +++ b/process/main.py @@ -94,7 +94,7 @@ ) from process.log import logging_model_handler, show_errors from process.pfcoil import PFCoil -from process.physics import DetailedPhysics, Physics +from process.physics import DetailedPhysics, Physics, PlasmaBeta from process.plasma_geometry import PlasmaGeom from process.plasma_profiles import PlasmaProfile from process.power import Power @@ -682,8 +682,11 @@ def __init__(self): neutral_beam=NeutralBeam(plasma_profile=self.plasma_profile), electron_bernstein=ElectronBernstein(plasma_profile=self.plasma_profile), ) + self.plasma_beta = PlasmaBeta() self.physics = Physics( - plasma_profile=self.plasma_profile, current_drive=self.current_drive + plasma_profile=self.plasma_profile, + current_drive=self.current_drive, + plasma_beta=self.plasma_beta, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/physics.py b/process/physics.py index 288be60f3..f622a6285 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1566,11 +1566,12 @@ def _trapped_particle_fraction_sauter( class Physics: - def __init__(self, plasma_profile, current_drive): + def __init__(self, plasma_profile, current_drive, plasma_beta): self.outfile = constants.NOUT self.mfile = constants.MFILE self.plasma_profile = plasma_profile self.current_drive = current_drive + self.beta = plasma_beta def physics(self): """ @@ -9034,6 +9035,14 @@ def reinke_tsep(b_plasma_toroidal_on_axis, flh, qstar, rmajor, eps, fgw, kappa, ) +class PlasmaBeta: + """Class to hold plasma beta calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" From 0f834861c4d428e5289021800e33d75b2704ac10 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 14:04:14 +0000 Subject: [PATCH 02/23] Refactor beta normalization calculations to use PlasmaBeta class methods --- process/physics.py | 298 +++++++++++++++++++++++---------------------- 1 file changed, 151 insertions(+), 147 deletions(-) diff --git a/process/physics.py b/process/physics.py index f622a6285..7773d8d58 100644 --- a/process/physics.py +++ b/process/physics.py @@ -2664,30 +2664,34 @@ def physics(self): # Define beta_norm_max calculations - physics_variables.beta_norm_max_wesson = self.calculate_beta_norm_max_wesson( - ind_plasma_internal_norm=physics_variables.ind_plasma_internal_norm + physics_variables.beta_norm_max_wesson = ( + self.beta.calculate_beta_norm_max_wesson( + ind_plasma_internal_norm=physics_variables.ind_plasma_internal_norm + ) ) # Original scaling law physics_variables.beta_norm_max_original_scaling = ( - self.calculate_beta_norm_max_original(eps=physics_variables.eps) + self.beta.calculate_beta_norm_max_original(eps=physics_variables.eps) ) # J. Menard scaling law - physics_variables.beta_norm_max_menard = self.calculate_beta_norm_max_menard( - eps=physics_variables.eps + physics_variables.beta_norm_max_menard = ( + self.beta.calculate_beta_norm_max_menard(eps=physics_variables.eps) ) # E. Tholerus scaling law - physics_variables.beta_norm_max_thloreus = self.calculate_beta_norm_max_thloreus( - c_beta=physics_variables.c_beta, - pres_plasma_on_axis=physics_variables.pres_plasma_thermal_on_axis, - pres_plasma_vol_avg=physics_variables.pres_plasma_thermal_vol_avg, + physics_variables.beta_norm_max_thloreus = ( + self.beta.calculate_beta_norm_max_thloreus( + c_beta=physics_variables.c_beta, + pres_plasma_on_axis=physics_variables.pres_plasma_thermal_on_axis, + pres_plasma_vol_avg=physics_variables.pres_plasma_thermal_vol_avg, + ) ) # R. D. Stambaugh scaling law physics_variables.beta_norm_max_stambaugh = ( - self.calculate_beta_norm_max_stambaugh( + self.beta.calculate_beta_norm_max_stambaugh( f_c_plasma_bootstrap=current_drive_variables.f_c_plasma_bootstrap, kappa=physics_variables.kappa, aspect=physics_variables.aspect, @@ -2953,143 +2957,6 @@ def calculate_internal_inductance_wesson(alphaj: float) -> float: """ return np.log(1.65 + 0.89 * alphaj) - @staticmethod - def calculate_beta_norm_max_wesson(ind_plasma_internal_norm: float) -> float: - """ - Calculate the Wesson normalsied beta upper limit. - - :param ind_plasma_internal_norm: Plasma normalised internal inductance - :type ind_plasma_internal_norm: float - - :return: The Wesson normalised beta upper limit. - :rtype: float - - :Notes: - - It is recommended to use this method with the other Wesson relations for normalsied internal - inductance and current profile index. - - This fit is derived from the DIII-D database for β_N >= 2.5 - - :References: - - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, - International Series of Monographs on Physics, Volume 149. - - - T. T. S et al., “Profile Optimization and High Beta Discharges and Stability of High Elongation Plasmas in the DIII-D Tokamak,” - Osti.gov, Oct. 1990. https://www.osti.gov/biblio/6194284 (accessed Dec. 19, 2024). - """ - return 4 * ind_plasma_internal_norm - - @staticmethod - def calculate_beta_norm_max_original(eps: float) -> float: - """ - Calculate the original scaling law normalsied beta upper limit. - - :param eps: Plasma normalised internal inductance - :type eps: float - - :return: The original scaling law normalised beta upper limit. - :rtype: float - - :Notes: - - :References: - - """ - return 2.7 * (1.0 + 5.0 * eps**3.5) - - @staticmethod - def calculate_beta_norm_max_menard(eps: float) -> float: - """ - Calculate the Menard normalsied beta upper limit. - - :param eps: Plasma normalised internal inductance - :type eps: float - - :return: The Menard normalised beta upper limit. - :rtype: float - - :Notes: - - Found as a reasonable fit to the computed no wall limit at f_BS ≈ 50% - - Uses maximum κ data from NSTX at A = 1.45, A = 1.75. Along with record - β_T data from DIII-D at A = 2.9 and high κ. - - :References: - - # J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” - Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, - doi: https://doi.org/10.1088/0029-5515/56/10/106023. - - """ - return 3.12 + 3.5 * eps**1.7 - - @staticmethod - def calculate_beta_norm_max_thloreus( - c_beta: float, pres_plasma_on_axis: float, pres_plasma_vol_avg: float - ) -> float: - """ - Calculate the E. Tholerus normalized beta upper limit. - - :param c_beta: Pressure peaking factor coefficient. - :type c_beta: float - :param pres_plasma_on_axis: Central plasma pressure (Pa). - :type pres_plasma_on_axis: float - :param pres_plasma_vol_avg: Volume-averaged plasma pressure (Pa). - :type pres_plasma_vol_avg: float - - :return: The E. Tholerus normalized beta upper limit. - :rtype: float - - :Notes: - - This method calculates the normalized beta upper limit based on the pressure peaking factor (Fp), - which is defined as the ratio of the peak pressure to the average pressure. - - The formula is derived from operational space studies of flat-top plasma in the STEP power plant. - - :References: - - E. Tholerus et al., “Flat-top plasma operational space of the STEP power plant,” - Nuclear Fusion, Aug. 2024, doi: https://doi.org/10.1088/1741-4326/ad6ea2. - """ - return 3.7 + ( - (c_beta / (pres_plasma_on_axis / pres_plasma_vol_avg)) - * (12.5 - 3.5 * (pres_plasma_on_axis / pres_plasma_vol_avg)) - ) - - @staticmethod - def calculate_beta_norm_max_stambaugh( - f_c_plasma_bootstrap: float, - kappa: float, - aspect: float, - ) -> float: - """ - Calculate the Stambaugh normalized beta upper limit. - - :param f_c_plasma_bootstrap: Bootstrap current fraction. - :type f_c_plasma_bootstrap: float - :param kappa: Plasma separatrix elongation. - :type kappa: float - :param aspect: Plasma aspect ratio. - :type aspect: float - - :return: The Stambaugh normalized beta upper limit. - :rtype: float - - :Notes: - - This method calculates the normalized beta upper limit based on the Stambaugh scaling. - - The formula is derived from empirical fits to high-performance, steady-state tokamak equilibria. - - :References: - - R. D. Stambaugh et al., “Fusion Nuclear Science Facility Candidates,” - Fusion Science and Technology, vol. 59, no. 2, pp. 279-307, Feb. 2011, - doi: https://doi.org/10.13182/fst59-279. - - - Y. R. Lin-Liu and R. D. Stambaugh, “Optimum equilibria for high performance, steady state tokamaks,” - Nuclear Fusion, vol. 44, no. 4, pp. 548-554, Mar. 2004, - doi: https://doi.org/10.1088/0029-5515/44/4/009. - """ - return ( - f_c_plasma_bootstrap - * 10 - * (-0.7748 + (1.2869 * kappa) - (0.2921 * kappa**2) + (0.0197 * kappa**3)) - / (aspect**0.5523 * np.tanh((1.8524 + (0.2319 * kappa)) / aspect**0.6163)) - ) - @staticmethod def calculate_internal_inductance_menard(kappa: float) -> float: """ @@ -9042,6 +8909,143 @@ def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE + @staticmethod + def calculate_beta_norm_max_wesson(ind_plasma_internal_norm: float) -> float: + """ + Calculate the Wesson normalsied beta upper limit. + + :param ind_plasma_internal_norm: Plasma normalised internal inductance + :type ind_plasma_internal_norm: float + + :return: The Wesson normalised beta upper limit. + :rtype: float + + :Notes: + - It is recommended to use this method with the other Wesson relations for normalsied internal + inductance and current profile index. + - This fit is derived from the DIII-D database for β_N >= 2.5 + + :References: + - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, + International Series of Monographs on Physics, Volume 149. + + - T. T. S et al., “Profile Optimization and High Beta Discharges and Stability of High Elongation Plasmas in the DIII-D Tokamak,” + Osti.gov, Oct. 1990. https://www.osti.gov/biblio/6194284 (accessed Dec. 19, 2024). + """ + return 4 * ind_plasma_internal_norm + + @staticmethod + def calculate_beta_norm_max_original(eps: float) -> float: + """ + Calculate the original scaling law normalsied beta upper limit. + + :param eps: Plasma normalised internal inductance + :type eps: float + + :return: The original scaling law normalised beta upper limit. + :rtype: float + + :Notes: + + :References: + + """ + return 2.7 * (1.0 + 5.0 * eps**3.5) + + @staticmethod + def calculate_beta_norm_max_menard(eps: float) -> float: + """ + Calculate the Menard normalsied beta upper limit. + + :param eps: Plasma normalised internal inductance + :type eps: float + + :return: The Menard normalised beta upper limit. + :rtype: float + + :Notes: + - Found as a reasonable fit to the computed no wall limit at f_BS ≈ 50% + - Uses maximum κ data from NSTX at A = 1.45, A = 1.75. Along with record + β_T data from DIII-D at A = 2.9 and high κ. + + :References: + - # J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” + Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, + doi: https://doi.org/10.1088/0029-5515/56/10/106023. + + """ + return 3.12 + 3.5 * eps**1.7 + + @staticmethod + def calculate_beta_norm_max_thloreus( + c_beta: float, pres_plasma_on_axis: float, pres_plasma_vol_avg: float + ) -> float: + """ + Calculate the E. Tholerus normalized beta upper limit. + + :param c_beta: Pressure peaking factor coefficient. + :type c_beta: float + :param pres_plasma_on_axis: Central plasma pressure (Pa). + :type pres_plasma_on_axis: float + :param pres_plasma_vol_avg: Volume-averaged plasma pressure (Pa). + :type pres_plasma_vol_avg: float + + :return: The E. Tholerus normalized beta upper limit. + :rtype: float + + :Notes: + - This method calculates the normalized beta upper limit based on the pressure peaking factor (Fp), + which is defined as the ratio of the peak pressure to the average pressure. + - The formula is derived from operational space studies of flat-top plasma in the STEP power plant. + + :References: + - E. Tholerus et al., “Flat-top plasma operational space of the STEP power plant,” + Nuclear Fusion, Aug. 2024, doi: https://doi.org/10.1088/1741-4326/ad6ea2. + """ + return 3.7 + ( + (c_beta / (pres_plasma_on_axis / pres_plasma_vol_avg)) + * (12.5 - 3.5 * (pres_plasma_on_axis / pres_plasma_vol_avg)) + ) + + @staticmethod + def calculate_beta_norm_max_stambaugh( + f_c_plasma_bootstrap: float, + kappa: float, + aspect: float, + ) -> float: + """ + Calculate the Stambaugh normalized beta upper limit. + + :param f_c_plasma_bootstrap: Bootstrap current fraction. + :type f_c_plasma_bootstrap: float + :param kappa: Plasma separatrix elongation. + :type kappa: float + :param aspect: Plasma aspect ratio. + :type aspect: float + + :return: The Stambaugh normalized beta upper limit. + :rtype: float + + :Notes: + - This method calculates the normalized beta upper limit based on the Stambaugh scaling. + - The formula is derived from empirical fits to high-performance, steady-state tokamak equilibria. + + :References: + - R. D. Stambaugh et al., “Fusion Nuclear Science Facility Candidates,” + Fusion Science and Technology, vol. 59, no. 2, pp. 279-307, Feb. 2011, + doi: https://doi.org/10.13182/fst59-279. + + - Y. R. Lin-Liu and R. D. Stambaugh, “Optimum equilibria for high performance, steady state tokamaks,” + Nuclear Fusion, vol. 44, no. 4, pp. 548-554, Mar. 2004, + doi: https://doi.org/10.1088/0029-5515/44/4/009. + """ + return ( + f_c_plasma_bootstrap + * 10 + * (-0.7748 + (1.2869 * kappa) - (0.2921 * kappa**2) + (0.0197 * kappa**3)) + / (aspect**0.5523 * np.tanh((1.8524 + (0.2319 * kappa)) / aspect**0.6163)) + ) + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" From dac541f317f476f802bdd73617c48394af146de4 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 14:15:14 +0000 Subject: [PATCH 03/23] Add static method to PlasmaBeta for calculating normalized beta --- process/physics.py | 24 ++++++++++++++++++++++++ tests/unit/test_physics.py | 25 ++++++++++++++++++++----- 2 files changed, 44 insertions(+), 5 deletions(-) diff --git a/process/physics.py b/process/physics.py index 7773d8d58..3bf5062b6 100644 --- a/process/physics.py +++ b/process/physics.py @@ -9046,6 +9046,30 @@ def calculate_beta_norm_max_stambaugh( / (aspect**0.5523 * np.tanh((1.8524 + (0.2319 * kappa)) / aspect**0.6163)) ) + @staticmethod + def calculate_normalised_beta( + beta: float, rminor: float, c_plasma: float, b_field: float + ) -> float: + """Calculate normalised beta (β_N). + + :param beta: Plasma beta (fraction). + :type beta: float + :param rminor: Plasma minor radius (m). + :type rminor: float + :param c_plasma: Plasma current (A). + :type c_plasma: float + :param b_field: Magnetic field (T). + :type b_field: float + :return: Normalised beta. + :rtype: float + + :Notes: + - 1.0e8 is a conversion factor to get beta_N in standard units, as plasma current is normally in MA and + beta is in percentage instead of fraction. + """ + + return 1.0e8 * (beta * rminor * b_field) / c_plasma + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 3caf87498..e1163391f 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -23,6 +23,7 @@ from process.physics import ( DetailedPhysics, Physics, + PlasmaBeta, calculate_beta_limit, calculate_current_coefficient_hastie, calculate_plasma_current_peng, @@ -55,6 +56,7 @@ def physics(): electron_bernstein=ElectronBernstein(plasma_profile=PlasmaProfile()), lower_hybrid=LowerHybrid(plasma_profile=PlasmaProfile()), ), + PlasmaBeta(), ) @@ -3421,21 +3423,21 @@ def test_calculate_internal_inductance_menard(): def test_calculate_beta_norm_max_wesson(): """Test calculate_beta_norm_max_wesson().""" ind_plasma_internal_norm = 1.5 - result = Physics.calculate_beta_norm_max_wesson(ind_plasma_internal_norm) + result = PlasmaBeta.calculate_beta_norm_max_wesson(ind_plasma_internal_norm) assert result == pytest.approx(6.0, abs=0.001) def test_calculate_beta_norm_max_original(): """Test calculate_beta_norm_max_original()""" eps = 0.5 - result = Physics.calculate_beta_norm_max_original(eps) + result = PlasmaBeta.calculate_beta_norm_max_original(eps) assert result == pytest.approx(3.8932426932522994, abs=0.00001) def test_calculate_beta_norm_max_menard(): """Test calculate_beta_norm_max_menard().""" eps = 0.5 - result = Physics.calculate_beta_norm_max_menard(eps) + result = PlasmaBeta.calculate_beta_norm_max_menard(eps) assert result == pytest.approx(4.197251361676802, abs=0.000001) @@ -3444,7 +3446,7 @@ def test_calculate_beta_norm_max_thloreus(): c_beta = 0.5 pres_plasma_on_axis = 2.0 pres_plasma_vol_avg = 1.0 - result = Physics.calculate_beta_norm_max_thloreus( + result = PlasmaBeta.calculate_beta_norm_max_thloreus( c_beta, pres_plasma_on_axis, pres_plasma_vol_avg ) assert result == pytest.approx(5.075, abs=0.00001) @@ -3455,7 +3457,7 @@ def test_calculate_beta_norm_max_stambaugh(): f_c_plasma_bootstrap = 0.7 kappa = 2.0 aspect = 2.5 - result = Physics.calculate_beta_norm_max_stambaugh( + result = PlasmaBeta.calculate_beta_norm_max_stambaugh( f_c_plasma_bootstrap, kappa, aspect ) assert result == pytest.approx(3.840954484207041, abs=0.00001) @@ -3469,6 +3471,19 @@ def test_calculate_internal_inductance_iter_3(): assert result == pytest.approx(0.9078959099585583, abs=0.00001) +def test_calculate_normalised_beta(): + """Test calculate_normalised_beta()""" + beta = 0.05 + rminor = 2.0 + c_plasma = 15.0e6 + b_field = 5.0 + + result = PlasmaBeta.calculate_normalised_beta( + beta=beta, rminor=rminor, c_plasma=c_plasma, b_field=b_field + ) + assert result == pytest.approx(3.3333333333333335, abs=1e-6) + + @pytest.mark.parametrize( "b_field, m_particle, z_particle, expected", [ From 4e188d92e67e7b5e73abaa6d08fbf101ecfb6ffe Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 14:42:31 +0000 Subject: [PATCH 04/23] Refactor beta information output to PlasmaBeta class for improved organization --- process/physics.py | 499 +++++++++++++++++---------------- tests/unit/test_stellarator.py | 3 +- 2 files changed, 255 insertions(+), 247 deletions(-) diff --git a/process/physics.py b/process/physics.py index 3bf5062b6..0df234eab 100644 --- a/process/physics.py +++ b/process/physics.py @@ -4381,253 +4381,9 @@ def outplas(self): po.oblnkl(self.outfile) po.ostars(self.outfile, 110) po.oblnkl(self.outfile) - po.osubhd(self.outfile, "Beta Information :") - if physics_variables.i_beta_component == 0: - po.ovarrf( - self.outfile, - "Upper limit on total beta", - "(beta_vol_avg_max)", - physics_variables.beta_vol_avg_max, - "OP ", - ) - elif physics_variables.i_beta_component == 1: - po.ovarrf( - self.outfile, - "Upper limit on thermal beta", - "(beta_vol_avg_max)", - physics_variables.beta_vol_avg_max, - "OP ", - ) - else: - po.ovarrf( - self.outfile, - "Upper limit on thermal + NB beta", - "(beta_vol_avg_max)", - physics_variables.beta_vol_avg_max, - "OP ", - ) - - po.ovarre( - self.outfile, - "Total plasma beta", - "(beta_total_vol_avg)", - physics_variables.beta_total_vol_avg, - ) - if physics_variables.i_beta_component == 0: - po.ovarrf( - self.outfile, - "Lower limit on total beta", - "(beta_vol_avg_min)", - physics_variables.beta_vol_avg_min, - "IP", - ) - elif physics_variables.i_beta_component == 1: - po.ovarrf( - self.outfile, - "Lower limit on thermal beta", - "(beta_vol_avg_min)", - physics_variables.beta_vol_avg_min, - "IP", - ) - else: - po.ovarrf( - self.outfile, - "Lower limit on thermal + NB beta", - "(beta_vol_avg_min)", - physics_variables.beta_vol_avg_min, - "IP", - ) - po.ovarre( - self.outfile, - "Upper limit on poloidal beta", - "(beta_poloidal_max)", - constraint_variables.beta_poloidal_max, - "IP", - ) - po.ovarre( - self.outfile, - "Total poloidal beta", - "(beta_poloidal_vol_avg)", - physics_variables.beta_poloidal_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Volume averaged toroidal beta", - "(beta_toroidal_vol_avg)", - physics_variables.beta_toroidal_vol_avg, - "OP ", - ) - for i in range(len(physics_variables.beta_thermal_toroidal_profile)): - po.ovarre( - self.mfile, - f"Beta toroidal profile at point {i}", - f"beta_thermal_toroidal_profile{i}", - physics_variables.beta_thermal_toroidal_profile[i], - ) - - po.ovarre( - self.outfile, - "Fast alpha beta", - "(beta_fast_alpha)", - physics_variables.beta_fast_alpha, - "OP ", - ) - po.ovarre( - self.outfile, - "Neutral Beam ion beta", - "(beta_beam)", - physics_variables.beta_beam, - "OP ", - ) - po.ovarre( - self.outfile, - "Ratio of fast alpha and beam beta to thermal beta", - "(f_beta_alpha_beam_thermal)", - physics_variables.f_beta_alpha_beam_thermal, - "OP ", - ) - - po.ovarre( - self.outfile, - "Volume averaged thermal beta", - "(beta_thermal_vol_avg)", - physics_variables.beta_thermal_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Thermal poloidal beta", - "(beta_thermal_poloidal_vol_avg)", - physics_variables.beta_thermal_poloidal_vol_avg, - "OP ", - ) - po.ovarre( - self.outfile, - "Thermal toroidal beta", - "(beta_thermal_toroidal_vol_avg)", - physics_variables.beta_thermal_toroidal_vol_avg, - "OP ", - ) - - po.ovarrf( - self.outfile, - "Poloidal beta and inverse aspect ratio", - "(beta_poloidal_eps)", - physics_variables.beta_poloidal_eps, - "OP ", - ) - po.ovarrf( - self.outfile, - "Poloidal beta and inverse aspect ratio upper limit", - "(beta_poloidal_eps_max)", - physics_variables.beta_poloidal_eps_max, - ) - po.osubhd(self.outfile, "Normalised Beta Information :") - if stellarator_variables.istell == 0: - if physics_variables.i_beta_norm_max != 0: - po.ovarrf( - self.outfile, - "Beta g coefficient", - "(beta_norm_max)", - physics_variables.beta_norm_max, - "OP ", - ) - else: - po.ovarrf( - self.outfile, - "Beta g coefficient", - "(beta_norm_max)", - physics_variables.beta_norm_max, - ) - po.ovarrf( - self.outfile, - "Normalised total beta", - "(beta_norm_total)", - physics_variables.beta_norm_total, - "OP ", - ) - po.ovarrf( - self.outfile, - "Normalised thermal beta", - "(beta_norm_thermal) ", - physics_variables.beta_norm_thermal, - "OP ", - ) - - po.ovarrf( - self.outfile, - "Normalised toroidal beta", - "(beta_norm_toroidal) ", - physics_variables.beta_norm_toroidal, - "OP ", - ) - - po.ovarrf( - self.outfile, - "Normalised poloidal beta", - "(beta_norm_poloidal) ", - physics_variables.beta_norm_poloidal, - "OP ", - ) - - po.osubhd(self.outfile, "Maximum normalised beta scalings :") - po.ovarrf( - self.outfile, - "J. Wesson normalised beta upper limit", - "(beta_norm_max_wesson) ", - physics_variables.beta_norm_max_wesson, - "OP ", - ) - po.ovarrf( - self.outfile, - "Original normalsied beta upper limit", - "(beta_norm_max_original_scaling) ", - physics_variables.beta_norm_max_original_scaling, - "OP ", - ) - po.ovarrf( - self.outfile, - "J. Menard normalised beta upper limit", - "(beta_norm_max_menard) ", - physics_variables.beta_norm_max_menard, - "OP ", - ) - po.ovarrf( - self.outfile, - "E. Thloreus normalised beta upper limit", - "(beta_norm_max_thloreus) ", - physics_variables.beta_norm_max_thloreus, - "OP ", - ) - po.ovarrf( - self.outfile, - "R. Stambaugh normalised beta upper limit", - "(beta_norm_max_stambaugh) ", - physics_variables.beta_norm_max_stambaugh, - "OP ", - ) - po.osubhd(self.outfile, "Plasma energies derived from beta :") - po.ovarre( - self.outfile, - "Plasma thermal energy derived from thermal beta (J)", - "(e_plasma_beta_thermal) ", - physics_variables.e_plasma_beta_thermal, - "OP", - ) - - po.ovarre( - self.outfile, - "Plasma thermal energy derived from the total beta (J)", - "(e_plasma_beta)", - physics_variables.e_plasma_beta, - "OP", - ) - - po.oblnkl(self.outfile) - po.ostars(self.outfile, 110) - po.oblnkl(self.outfile) + # Output beta information + self.beta.output_beta_information() po.osubhd(self.outfile, "Temperature and Density (volume averaged) :") po.ovarrf( @@ -9070,6 +8826,257 @@ def calculate_normalised_beta( return 1.0e8 * (beta * rminor * b_field) / c_plasma + def output_beta_information(self): + """Output beta information to file.""" + + po.osubhd(self.outfile, "Beta Information :") + if physics_variables.i_beta_component == 0: + po.ovarrf( + self.outfile, + "Upper limit on total beta", + "(beta_vol_avg_max)", + physics_variables.beta_vol_avg_max, + "OP ", + ) + elif physics_variables.i_beta_component == 1: + po.ovarrf( + self.outfile, + "Upper limit on thermal beta", + "(beta_vol_avg_max)", + physics_variables.beta_vol_avg_max, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Upper limit on thermal + NB beta", + "(beta_vol_avg_max)", + physics_variables.beta_vol_avg_max, + "OP ", + ) + + po.ovarre( + self.outfile, + "Total plasma beta", + "(beta_total_vol_avg)", + physics_variables.beta_total_vol_avg, + ) + if physics_variables.i_beta_component == 0: + po.ovarrf( + self.outfile, + "Lower limit on total beta", + "(beta_vol_avg_min)", + physics_variables.beta_vol_avg_min, + "IP", + ) + elif physics_variables.i_beta_component == 1: + po.ovarrf( + self.outfile, + "Lower limit on thermal beta", + "(beta_vol_avg_min)", + physics_variables.beta_vol_avg_min, + "IP", + ) + else: + po.ovarrf( + self.outfile, + "Lower limit on thermal + NB beta", + "(beta_vol_avg_min)", + physics_variables.beta_vol_avg_min, + "IP", + ) + po.ovarre( + self.outfile, + "Upper limit on poloidal beta", + "(beta_poloidal_max)", + constraint_variables.beta_poloidal_max, + "IP", + ) + po.ovarre( + self.outfile, + "Total poloidal beta", + "(beta_poloidal_vol_avg)", + physics_variables.beta_poloidal_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Volume averaged toroidal beta", + "(beta_toroidal_vol_avg)", + physics_variables.beta_toroidal_vol_avg, + "OP ", + ) + for i in range(len(physics_variables.beta_thermal_toroidal_profile)): + po.ovarre( + self.mfile, + f"Beta toroidal profile at point {i}", + f"beta_thermal_toroidal_profile{i}", + physics_variables.beta_thermal_toroidal_profile[i], + ) + + po.ovarre( + self.outfile, + "Fast alpha beta", + "(beta_fast_alpha)", + physics_variables.beta_fast_alpha, + "OP ", + ) + po.ovarre( + self.outfile, + "Neutral Beam ion beta", + "(beta_beam)", + physics_variables.beta_beam, + "OP ", + ) + po.ovarre( + self.outfile, + "Ratio of fast alpha and beam beta to thermal beta", + "(f_beta_alpha_beam_thermal)", + physics_variables.f_beta_alpha_beam_thermal, + "OP ", + ) + + po.ovarre( + self.outfile, + "Volume averaged thermal beta", + "(beta_thermal_vol_avg)", + physics_variables.beta_thermal_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Thermal poloidal beta", + "(beta_thermal_poloidal_vol_avg)", + physics_variables.beta_thermal_poloidal_vol_avg, + "OP ", + ) + po.ovarre( + self.outfile, + "Thermal toroidal beta", + "(beta_thermal_toroidal_vol_avg)", + physics_variables.beta_thermal_toroidal_vol_avg, + "OP ", + ) + + po.ovarrf( + self.outfile, + "Poloidal beta and inverse aspect ratio", + "(beta_poloidal_eps)", + physics_variables.beta_poloidal_eps, + "OP ", + ) + po.ovarrf( + self.outfile, + "Poloidal beta and inverse aspect ratio upper limit", + "(beta_poloidal_eps_max)", + physics_variables.beta_poloidal_eps_max, + ) + po.osubhd(self.outfile, "Normalised Beta Information :") + if stellarator_variables.istell == 0: + if physics_variables.i_beta_norm_max != 0: + po.ovarrf( + self.outfile, + "Beta g coefficient", + "(beta_norm_max)", + physics_variables.beta_norm_max, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Beta g coefficient", + "(beta_norm_max)", + physics_variables.beta_norm_max, + ) + po.ovarrf( + self.outfile, + "Normalised total beta", + "(beta_norm_total)", + physics_variables.beta_norm_total, + "OP ", + ) + po.ovarrf( + self.outfile, + "Normalised thermal beta", + "(beta_norm_thermal) ", + physics_variables.beta_norm_thermal, + "OP ", + ) + + po.ovarrf( + self.outfile, + "Normalised toroidal beta", + "(beta_norm_toroidal) ", + physics_variables.beta_norm_toroidal, + "OP ", + ) + + po.ovarrf( + self.outfile, + "Normalised poloidal beta", + "(beta_norm_poloidal) ", + physics_variables.beta_norm_poloidal, + "OP ", + ) + + po.osubhd(self.outfile, "Maximum normalised beta scalings :") + po.ovarrf( + self.outfile, + "J. Wesson normalised beta upper limit", + "(beta_norm_max_wesson) ", + physics_variables.beta_norm_max_wesson, + "OP ", + ) + po.ovarrf( + self.outfile, + "Original normalsied beta upper limit", + "(beta_norm_max_original_scaling) ", + physics_variables.beta_norm_max_original_scaling, + "OP ", + ) + po.ovarrf( + self.outfile, + "J. Menard normalised beta upper limit", + "(beta_norm_max_menard) ", + physics_variables.beta_norm_max_menard, + "OP ", + ) + po.ovarrf( + self.outfile, + "E. Thloreus normalised beta upper limit", + "(beta_norm_max_thloreus) ", + physics_variables.beta_norm_max_thloreus, + "OP ", + ) + po.ovarrf( + self.outfile, + "R. Stambaugh normalised beta upper limit", + "(beta_norm_max_stambaugh) ", + physics_variables.beta_norm_max_stambaugh, + "OP ", + ) + + po.osubhd(self.outfile, "Plasma energies derived from beta :") + po.ovarre( + self.outfile, + "Plasma thermal energy derived from thermal beta (J)", + "(e_plasma_beta_thermal) ", + physics_variables.e_plasma_beta_thermal, + "OP", + ) + + po.ovarre( + self.outfile, + "Plasma thermal energy derived from the total beta (J)", + "(e_plasma_beta)", + physics_variables.e_plasma_beta, + "OP", + ) + + po.oblnkl(self.outfile) + po.ostars(self.outfile, 110) + po.oblnkl(self.outfile) + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index db5d4df06..44c7c77c9 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -28,7 +28,7 @@ ) from process.fw import Fw from process.hcpb import CCFE_HCPB -from process.physics import Physics +from process.physics import Physics, PlasmaBeta from process.plasma_profiles import PlasmaProfile from process.power import Power from process.stellarator import Neoclassics, Stellarator @@ -68,6 +68,7 @@ def stellarator(): LowerHybrid(plasma_profile=PlasmaProfile()), ElectronBernstein(plasma_profile=PlasmaProfile()), ), + PlasmaBeta(), ), Neoclassics(), ) From ce7c70069b203acbe9dc7be519f6ebee5115c471 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 14:49:27 +0000 Subject: [PATCH 05/23] Move total normalised beta calc to outside of calculate_plasma_current --- process/physics.py | 17 ++++++++--------- tests/unit/test_physics.py | 10 ---------- 2 files changed, 8 insertions(+), 19 deletions(-) diff --git a/process/physics.py b/process/physics.py index 0df234eab..d469425b3 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1644,6 +1644,14 @@ def physics(self): physics_variables.triang95, ) + # Normalised beta from Troyon beta limit + physics_variables.beta_norm_total = self.beta.calculate_normalised_beta( + beta=physics_variables.beta_total_vol_avg, + rminor=physics_variables.rminor, + c_plasma=physics_variables.plasma_current, + b_field=physics_variables.b_plasma_toroidal_on_axis, + ) + # ----------------------------------------------------- # Plasma Current Profile # ----------------------------------------------------- @@ -3808,15 +3816,6 @@ def calculate_plasma_current( * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) ) - # Normalised beta from Troyon beta limit - physics_variables.beta_norm_total = ( - 1.0e8 - * physics_variables.beta_total_vol_avg - * rminor - * b_plasma_toroidal_on_axis - / plasma_current - ) - # Calculate the poloidal field generated by the plasma current b_plasma_poloidal_average = calculate_poloidal_field( i_plasma_current, diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index e1163391f..6a2bf6b24 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -1215,12 +1215,6 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): :type monkeypatch: _pytest.monkeypatch.monkeypatch """ - monkeypatch.setattr( - physics_variables, - "beta_norm_total", - plasmacurrentparam.beta_norm_total, - ) - monkeypatch.setattr(physics_variables, "beta_total_vol_avg", plasmacurrentparam.beta) b_plasma_poloidal_average, qstar, plasma_current = physics.calculate_plasma_current( @@ -1240,10 +1234,6 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): triang95=plasmacurrentparam.triang95, ) - assert physics_variables.beta_norm_total == pytest.approx( - plasmacurrentparam.expected_normalised_total_beta - ) - assert b_plasma_poloidal_average == pytest.approx(plasmacurrentparam.expected_bp) assert qstar == pytest.approx(plasmacurrentparam.expected_qstar) From 303709003b218caea8760991284a0208e74cfbb0 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 15:17:54 +0000 Subject: [PATCH 06/23] Add static method to PlasmaBeta for calculating plasma energy from beta --- process/physics.py | 21 +++++++++++++++++++++ tests/unit/test_physics.py | 11 +++++++++++ 2 files changed, 32 insertions(+) diff --git a/process/physics.py b/process/physics.py index d469425b3..a792945f9 100644 --- a/process/physics.py +++ b/process/physics.py @@ -8825,6 +8825,27 @@ def calculate_normalised_beta( return 1.0e8 * (beta * rminor * b_field) / c_plasma + @staticmethod + def calculate_plasma_energy_from_beta( + beta: float, b_field: float, vol_plasma: float + ) -> float: + """Calculate plasma thermal energy from beta. + + E_plasma = 1.5 * β * B² / (2 * μ_0) * V + + :param beta: Plasma beta (fraction). + :type beta: float + :param b_field: Magnetic field (T). + :type b_field: float + :param vol_plasma: Plasma volume (m³). + :type vol_plasma: float + :return: Plasma energy (J). + :rtype: float + + """ + + return (1.5e0 * beta * b_field**2) / (2.0e0 * constants.RMU0) * vol_plasma + def output_beta_information(self): """Output beta information to file.""" diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 6a2bf6b24..d2b6f126a 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -3474,6 +3474,17 @@ def test_calculate_normalised_beta(): assert result == pytest.approx(3.3333333333333335, abs=1e-6) +def test_calculate_plasma_energy_from_beta(): + """Test calculate_plasma_energy_from_beta()""" + beta = 0.02 + b_field = 5.3 + vol_plasma = 1000.0 + result = PlasmaBeta.calculate_plasma_energy_from_beta( + beta=beta, b_field=b_field, vol_plasma=vol_plasma + ) + assert result == pytest.approx(335299676.2083403, rel=0.01) + + @pytest.mark.parametrize( "b_field, m_particle, z_particle, expected", [ From 54fd07048e6b477df4734b8127cabec2b4406176 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 15:54:17 +0000 Subject: [PATCH 07/23] Refactor plasma energy calculations to use PlasmaBeta methods for improved clarity and maintainability --- process/physics.py | 22 +++++++++------------- 1 file changed, 9 insertions(+), 13 deletions(-) diff --git a/process/physics.py b/process/physics.py index a792945f9..f0362a4ef 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1881,22 +1881,18 @@ def physics(self): # Plasma thermal energy derived from the thermal beta physics_variables.e_plasma_beta_thermal = ( - 1.5e0 - * physics_variables.beta_thermal_vol_avg - * physics_variables.b_plasma_total - * physics_variables.b_plasma_total - / (2.0e0 * constants.RMU0) - * physics_variables.vol_plasma + self.beta.calculate_plasma_energy_from_beta( + beta=physics_variables.beta_thermal_vol_avg, + b_field=physics_variables.b_plasma_total, + vol_plasma=physics_variables.vol_plasma, + ) ) # Plasma thermal energy derived from the total beta - physics_variables.e_plasma_beta = ( - 1.5e0 - * physics_variables.beta_total_vol_avg - * physics_variables.b_plasma_total - * physics_variables.b_plasma_total - / (2.0e0 * constants.RMU0) - * physics_variables.vol_plasma + physics_variables.e_plasma_beta = self.beta.calculate_plasma_energy_from_beta( + beta=physics_variables.beta_total_vol_avg, + b_field=physics_variables.b_plasma_total, + vol_plasma=physics_variables.vol_plasma, ) # ======================================================= From 108a973b29e96c227df41cf718876cd55948fec6 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Tue, 20 Jan 2026 16:19:15 +0000 Subject: [PATCH 08/23] Refactor fast alpha beta calculation into PlasmaBeta class and remove redundant function from physics_functions.py --- process/main.py | 1 + process/physics.py | 139 ++++++++++++++++++++++++++++++--- process/physics_functions.py | 109 -------------------------- process/stellarator.py | 4 +- tests/unit/test_stellarator.py | 1 + 5 files changed, 133 insertions(+), 121 deletions(-) diff --git a/process/main.py b/process/main.py index 4d86437ee..5d3460782 100644 --- a/process/main.py +++ b/process/main.py @@ -703,6 +703,7 @@ def __init__(self): current_drive=self.current_drive, physics=self.physics, neoclassics=self.neoclassics, + plasma_beta=self.plasma_beta, ) self.dcll = DCLL(fw=self.fw) diff --git a/process/physics.py b/process/physics.py index f0362a4ef..284942a78 100644 --- a/process/physics.py +++ b/process/physics.py @@ -2354,17 +2354,17 @@ def physics(self): physics_variables.pden_plasma_alpha_mw, ) - physics_variables.beta_fast_alpha = physics_funcs.fast_alpha_beta( - physics_variables.b_plasma_poloidal_average, - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.nd_plasma_electrons_vol_avg, - physics_variables.nd_plasma_fuel_ions_vol_avg, - physics_variables.nd_plasma_ions_total_vol_avg, - physics_variables.temp_plasma_electron_density_weighted_kev, - physics_variables.temp_plasma_ion_density_weighted_kev, - physics_variables.pden_alpha_total_mw, - physics_variables.pden_plasma_alpha_mw, - physics_variables.i_beta_fast_alpha, + physics_variables.beta_fast_alpha = self.beta.fast_alpha_beta( + b_plasma_poloidal_average=physics_variables.b_plasma_poloidal_average, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + nd_plasma_electrons_vol_avg=physics_variables.nd_plasma_electrons_vol_avg, + nd_plasma_fuel_ions_vol_avg=physics_variables.nd_plasma_fuel_ions_vol_avg, + nd_plasma_ions_total_vol_avg=physics_variables.nd_plasma_ions_total_vol_avg, + temp_plasma_electron_density_weighted_kev=physics_variables.temp_plasma_electron_density_weighted_kev, + temp_plasma_ion_density_weighted_kev=physics_variables.temp_plasma_ion_density_weighted_kev, + pden_alpha_total_mw=physics_variables.pden_alpha_total_mw, + pden_plasma_alpha_mw=physics_variables.pden_plasma_alpha_mw, + i_beta_fast_alpha=physics_variables.i_beta_fast_alpha, ) # Nominal mean neutron wall load on entire first wall area including divertor and beam holes @@ -8842,6 +8842,123 @@ def calculate_plasma_energy_from_beta( return (1.5e0 * beta * b_field**2) / (2.0e0 * constants.RMU0) * vol_plasma + @staticmethod + def fast_alpha_beta( + b_plasma_poloidal_average: float, + b_plasma_toroidal_on_axis: float, + nd_plasma_electrons_vol_avg: float, + nd_plasma_fuel_ions_vol_avg: float, + nd_plasma_ions_total_vol_avg: float, + temp_plasma_electron_density_weighted_kev: float, + temp_plasma_ion_density_weighted_kev: float, + pden_alpha_total_mw: float, + pden_plasma_alpha_mw: float, + i_beta_fast_alpha: int, + ) -> float: + """ + Calculate the fast alpha beta component. + + This function computes the fast alpha beta contribution based on the provided plasma parameters. + + :param b_plasma_poloidal_average: Poloidal field (T). + :type b_plasma_poloidal_average: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param nd_plasma_electrons_vol_avg: Electron density (m^-3). + :type nd_plasma_electrons_vol_avg: float + :param nd_plasma_fuel_ions_vol_avg: Fuel ion density (m^-3). + :type nd_plasma_fuel_ions_vol_avg: float + :param nd_plasma_ions_total_vol_avg: Total ion density (m^-3). + :type nd_plasma_ions_total_vol_avg: float + :param temp_plasma_electron_density_weighted_kev: Density-weighted electron temperature (keV). + :type temp_plasma_electron_density_weighted_kev: float + :param temp_plasma_ion_density_weighted_kev: Density-weighted ion temperature (keV). + :type temp_plasma_ion_density_weighted_kev: float + :param pden_alpha_total_mw: Alpha power per unit volume, from beams and plasma (MW/m^3). + :type pden_alpha_total_mw: float + :param pden_plasma_alpha_mw: Alpha power per unit volume just from plasma (MW/m^3). + :type pden_plasma_alpha_mw: float + :param i_beta_fast_alpha: Switch for fast alpha pressure method. + :type i_beta_fast_alpha: int + + :return: Fast alpha beta component. + :rtype: float + + :Notes: + - For IPDG89 scaling applicability is Z_eff = 1.5, T_i/T_e = 1, = 5-20 keV + + + :References: + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', + https://inis.iaea.org/collection/NCLCollectionStore/_Public/21/068/21068960.pdf + + - Uckan, N. A., Tolliver, J. S., Houlberg, W. A., and Attenberger, S. E. + Influence of fast alpha diffusion and thermal alpha buildup on tokamak reactor performance. + United States: N. p., 1987. Web.https://www.osti.gov/servlets/purl/5611706 + + """ + + # Determine average fast alpha density + if physics_variables.f_plasma_fuel_deuterium < 1.0: + beta_thermal = ( + 2.0 + * constants.RMU0 + * constants.KILOELECTRON_VOLT + * ( + nd_plasma_electrons_vol_avg + * temp_plasma_electron_density_weighted_kev + + nd_plasma_ions_total_vol_avg * temp_plasma_ion_density_weighted_kev + ) + / (b_plasma_toroidal_on_axis**2 + b_plasma_poloidal_average**2) + ) + + # jlion: This "fact" model is heavily flawed for smaller temperatures! It is unphysical for a stellarator (high n low T) + # IPDG89 fast alpha scaling + if i_beta_fast_alpha == 0: + fact = min( + 0.3, + 0.29 + * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 + * ( + ( + temp_plasma_electron_density_weighted_kev + + temp_plasma_ion_density_weighted_kev + ) + / 20.0 + - 0.37 + ), + ) + + # Modified scaling, D J Ward + else: + fact = min( + 0.30, + 0.26 + * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 + * np.sqrt( + max( + 0.0, + ( + ( + temp_plasma_electron_density_weighted_kev + + temp_plasma_ion_density_weighted_kev + ) + / 20.0 + - 0.65 + ), + ) + ), + ) + + fact = max(fact, 0.0) + fact2 = pden_alpha_total_mw / pden_plasma_alpha_mw + beta_fast_alpha = beta_thermal * fact * fact2 + + else: # negligible alpha production, alpha_power_density = p_beam_alpha_mw = 0 + beta_fast_alpha = 0.0 + + return beta_fast_alpha + def output_beta_information(self): """Output beta information to file.""" diff --git a/process/physics_functions.py b/process/physics_functions.py index f0346fc5a..de5382776 100644 --- a/process/physics_functions.py +++ b/process/physics_functions.py @@ -4,8 +4,6 @@ import numpy as np import process.impurity_radiation as impurity -from process import constants -from process.data_structure import physics_variables from process.plasma_profiles import PlasmaProfile logger = logging.getLogger(__name__) @@ -219,110 +217,3 @@ def psync_albajar_fidone( # pden_plasma_sync_mw should be per unit volume; Albajar gives it as total return p_sync_mw / vol_plasma - - -def fast_alpha_beta( - b_plasma_poloidal_average: float, - b_plasma_toroidal_on_axis: float, - nd_plasma_electrons_vol_avg: float, - nd_plasma_fuel_ions_vol_avg: float, - nd_plasma_ions_total_vol_avg: float, - temp_plasma_electron_density_weighted_kev: float, - temp_plasma_ion_density_weighted_kev: float, - pden_alpha_total_mw: float, - pden_plasma_alpha_mw: float, - i_beta_fast_alpha: int, -) -> float: - """ - Calculate the fast alpha beta component. - - This function computes the fast alpha beta contribution based on the provided plasma parameters. - - Parameters: - b_plasma_poloidal_average (float): Poloidal field (T). - b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). - nd_plasma_electrons_vol_avg (float): Electron density (m^-3). - nd_plasma_fuel_ions_vol_avg (float): Fuel ion density (m^-3). - nd_plasma_ions_total_vol_avg (float): Total ion density (m^-3). - temp_plasma_electron_density_weighted_kev (float): Density-weighted electron temperature (keV). - temp_plasma_ion_density_weighted_kev (float): Density-weighted ion temperature (keV). - pden_alpha_total_mw (float): Alpha power per unit volume, from beams and plasma (MW/m^3). - pden_plasma_alpha_mw (float): Alpha power per unit volume just from plasma (MW/m^3). - i_beta_fast_alpha (int): Switch for fast alpha pressure method. - - Returns: - float: Fast alpha beta component. - - Notes: - - For IPDG89 scaling applicability is Z_eff = 1.5, T_i/T_e = 1, = 5-20 keV - - - References: - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', - https://inis.iaea.org/collection/NCLCollectionStore/_Public/21/068/21068960.pdf - - - Uckan, N. A., Tolliver, J. S., Houlberg, W. A., and Attenberger, S. E. - Influence of fast alpha diffusion and thermal alpha buildup on tokamak reactor performance. - United States: N. p., 1987. Web.https://www.osti.gov/servlets/purl/5611706 - - """ - - # Determine average fast alpha density - if physics_variables.f_plasma_fuel_deuterium < 1.0: - beta_thermal = ( - 2.0 - * constants.RMU0 - * constants.KILOELECTRON_VOLT - * ( - nd_plasma_electrons_vol_avg * temp_plasma_electron_density_weighted_kev - + nd_plasma_ions_total_vol_avg * temp_plasma_ion_density_weighted_kev - ) - / (b_plasma_toroidal_on_axis**2 + b_plasma_poloidal_average**2) - ) - - # jlion: This "fact" model is heavily flawed for smaller temperatures! It is unphysical for a stellarator (high n low T) - # IPDG89 fast alpha scaling - if i_beta_fast_alpha == 0: - fact = min( - 0.3, - 0.29 - * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 - * ( - ( - temp_plasma_electron_density_weighted_kev - + temp_plasma_ion_density_weighted_kev - ) - / 20.0 - - 0.37 - ), - ) - - # Modified scaling, D J Ward - else: - fact = min( - 0.30, - 0.26 - * (nd_plasma_fuel_ions_vol_avg / nd_plasma_electrons_vol_avg) ** 2 - * np.sqrt( - max( - 0.0, - ( - ( - temp_plasma_electron_density_weighted_kev - + temp_plasma_ion_density_weighted_kev - ) - / 20.0 - - 0.65 - ), - ) - ), - ) - - fact = max(fact, 0.0) - fact2 = pden_alpha_total_mw / pden_plasma_alpha_mw - beta_fast_alpha = beta_thermal * fact * fact2 - - else: # negligible alpha production, alpha_power_density = p_beam_alpha_mw = 0 - beta_fast_alpha = 0.0 - - return beta_fast_alpha diff --git a/process/stellarator.py b/process/stellarator.py index 34113d307..18533d959 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -66,6 +66,7 @@ def __init__( current_drive, physics, neoclassics, + plasma_beta, ) -> None: """Initialises the Stellarator model's variables @@ -102,6 +103,7 @@ def __init__( self.current_drive = current_drive self.physics = physics self.neoclassics = neoclassics + self.beta = plasma_beta def run(self, output: bool): """Routine to call the physics and engineering modules @@ -4382,7 +4384,7 @@ def stphys(self, output): physics_variables.pden_plasma_alpha_mw, ) - physics_variables.beta_fast_alpha = physics_funcs.fast_alpha_beta( + physics_variables.beta_fast_alpha = self.beta.fast_alpha_beta( physics_variables.b_plasma_poloidal_average, physics_variables.b_plasma_toroidal_on_axis, physics_variables.nd_plasma_electrons_vol_avg, diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 44c7c77c9..4d3320414 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -71,6 +71,7 @@ def stellarator(): PlasmaBeta(), ), Neoclassics(), + plasma_beta=PlasmaBeta(), ) From ec47444923c7598f6d6b3a2fdb50d47434ca3141 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 21 Jan 2026 08:55:58 +0000 Subject: [PATCH 09/23] Refactor beta limit calculation into PlasmaBeta class and remove redundant function from physics.py --- process/physics.py | 125 ++++++++++++++++++------------------- tests/unit/test_physics.py | 5 +- 2 files changed, 65 insertions(+), 65 deletions(-) diff --git a/process/physics.py b/process/physics.py index 284942a78..00cbf36c8 100644 --- a/process/physics.py +++ b/process/physics.py @@ -201,57 +201,6 @@ def calculate_volt_second_requirements( ) -@nb.jit(nopython=True, cache=True) -def calculate_beta_limit( - b_plasma_toroidal_on_axis: float, - beta_norm_max: float, - plasma_current: float, - rminor: float, -) -> float: - """ - Calculate the beta scaling limit. - - :param b_plasma_toroidal_on_axis: Toroidal B-field on plasma axis [T]. - :type b_plasma_toroidal_on_axis: float - :param beta_norm_max: Troyon-like g coefficient. - :type beta_norm_max: float - :param plasma_current: Plasma current [A]. - :type plasma_current: float - :param rminor: Plasma minor axis [m]. - :type rminor: float - :return: Beta limit as defined below. - :rtype: float - - This subroutine calculates the beta limit using the algorithm documented in AEA FUS 172. - The limit applies to beta defined with respect to the total B-field. - Switch i_beta_component determines which components of beta to include. - - Notes: - - If i_beta_component = 0, then the limit is applied to the total beta. - - If i_beta_component = 1, then the limit is applied to the thermal beta only. - - If i_beta_component = 2, then the limit is applied to the thermal + neutral beam beta components. - - If i_beta_component = 3, then the limit is applied to the toroidal beta. - - - The default value for the g coefficient is beta_norm_max = 3.5. - - References: - - F. Troyon et.al, “Beta limit in tokamaks. Experimental and computational status,” - Plasma Physics and Controlled Fusion, vol. 30, no. 11, pp. 1597-1609, Oct. 1988, - doi: https://doi.org/10.1088/0741-3335/30/11/019. - - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - - # Multiplied by 0.01 to convert from % to fraction - return ( - 0.01 - * beta_norm_max - * (plasma_current / 1.0e6) - / (rminor * b_plasma_toroidal_on_axis) - ) - - # ----------------------------------------------------- # Plasma Current & Poloidal Field Calculations # ----------------------------------------------------- @@ -2724,11 +2673,11 @@ def physics(self): ) # calculate_beta_limit() returns the beta_vol_avg_max for beta - physics_variables.beta_vol_avg_max = calculate_beta_limit( - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.beta_norm_max, - physics_variables.plasma_current, - physics_variables.rminor, + physics_variables.beta_vol_avg_max = self.beta.calculate_beta_limit_from_norm( + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + beta_norm_max=physics_variables.beta_norm_max, + plasma_current=physics_variables.plasma_current, + rminor=physics_variables.rminor, ) # ============================================================ @@ -8842,6 +8791,56 @@ def calculate_plasma_energy_from_beta( return (1.5e0 * beta * b_field**2) / (2.0e0 * constants.RMU0) * vol_plasma + @staticmethod + def calculate_beta_limit_from_norm( + b_plasma_toroidal_on_axis: float, + beta_norm_max: float, + plasma_current: float, + rminor: float, + ) -> float: + """ + Calculate the maximum allowed beta (β) from a given normalised (β_N). + + :param b_plasma_toroidal_on_axis: Toroidal B-field on plasma axis [T]. + :type b_plasma_toroidal_on_axis: float + :param beta_norm_max: Troyon-like g coefficient. + :type beta_norm_max: float + :param plasma_current: Plasma current [A]. + :type plasma_current: float + :param rminor: Plasma minor axis [m]. + :type rminor: float + :return: Beta limit as defined below. + :rtype: float + + This subroutine calculates the beta limit using the algorithm documented in AEA FUS 172. + The limit applies to beta defined with respect to the total B-field. + Switch i_beta_component determines which components of beta to include. + + Notes: + - If i_beta_component = 0, then the limit is applied to the total beta. + - If i_beta_component = 1, then the limit is applied to the thermal beta only. + - If i_beta_component = 2, then the limit is applied to the thermal + neutral beam beta components. + - If i_beta_component = 3, then the limit is applied to the toroidal beta. + + - The default value for the g coefficient is beta_norm_max = 3.5. + + References: + - F. Troyon et.al, “Beta limit in tokamaks. Experimental and computational status,” + Plasma Physics and Controlled Fusion, vol. 30, no. 11, pp. 1597-1609, Oct. 1988, + doi: https://doi.org/10.1088/0741-3335/30/11/019. + + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + + # Multiplied by 0.01 to convert from % to fraction + return ( + 0.01 + * beta_norm_max + * (plasma_current / 1.0e6) + / (rminor * b_plasma_toroidal_on_axis) + ) + @staticmethod def fast_alpha_beta( b_plasma_poloidal_average: float, @@ -8856,7 +8855,7 @@ def fast_alpha_beta( i_beta_fast_alpha: int, ) -> float: """ - Calculate the fast alpha beta component. + Calculate the fast alpha beta (β_fast_alpha) component. This function computes the fast alpha beta contribution based on the provided plasma parameters. @@ -8864,19 +8863,19 @@ def fast_alpha_beta( :type b_plasma_poloidal_average: float :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). :type b_plasma_toroidal_on_axis: float - :param nd_plasma_electrons_vol_avg: Electron density (m^-3). + :param nd_plasma_electrons_vol_avg: Electron density (m⁻³). :type nd_plasma_electrons_vol_avg: float - :param nd_plasma_fuel_ions_vol_avg: Fuel ion density (m^-3). + :param nd_plasma_fuel_ions_vol_avg: Fuel ion density (m⁻³). :type nd_plasma_fuel_ions_vol_avg: float - :param nd_plasma_ions_total_vol_avg: Total ion density (m^-3). + :param nd_plasma_ions_total_vol_avg: Total ion density (m⁻³). :type nd_plasma_ions_total_vol_avg: float :param temp_plasma_electron_density_weighted_kev: Density-weighted electron temperature (keV). :type temp_plasma_electron_density_weighted_kev: float :param temp_plasma_ion_density_weighted_kev: Density-weighted ion temperature (keV). :type temp_plasma_ion_density_weighted_kev: float - :param pden_alpha_total_mw: Alpha power per unit volume, from beams and plasma (MW/m^3). + :param pden_alpha_total_mw: Alpha power per unit volume, from beams and plasma (MW/m³). :type pden_alpha_total_mw: float - :param pden_plasma_alpha_mw: Alpha power per unit volume just from plasma (MW/m^3). + :param pden_plasma_alpha_mw: Alpha power per unit volume just from plasma (MW/m³). :type pden_plasma_alpha_mw: float :param i_beta_fast_alpha: Switch for fast alpha pressure method. :type i_beta_fast_alpha: int @@ -8885,7 +8884,7 @@ def fast_alpha_beta( :rtype: float :Notes: - - For IPDG89 scaling applicability is Z_eff = 1.5, T_i/T_e = 1, = 5-20 keV + - For IPDG89 scaling applicability is Z_eff = 1.5, T_i/T_e = 1, 〈T〉 = 5-20 keV :References: diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index d2b6f126a..e34f0d715 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -24,7 +24,6 @@ DetailedPhysics, Physics, PlasmaBeta, - calculate_beta_limit, calculate_current_coefficient_hastie, calculate_plasma_current_peng, calculate_poloidal_beta, @@ -1329,7 +1328,9 @@ def test_calculate_poloidal_field(arguments, expected): def test_calculate_beta_limit(): - assert calculate_beta_limit(12, 4.879, 18300000, 2.5) == pytest.approx(0.0297619) + assert PlasmaBeta.calculate_beta_limit_from_norm( + 12, 4.879, 18300000, 2.5 + ) == pytest.approx(0.0297619) def test_conhas(): From 8aac633b5761bcda656c74ee4743102e61fe38ea Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 21 Jan 2026 09:24:30 +0000 Subject: [PATCH 10/23] Refactor poloidal beta calculation to use PlasmaBeta class methods and remove redundant function --- process/physics.py | 44 ++++++++++++++++++++++---------------- tests/unit/test_physics.py | 3 +-- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/process/physics.py b/process/physics.py index 00cbf36c8..545a9230c 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1749,10 +1749,10 @@ def physics(self): ) # Calculate physics_variables.beta poloidal [-] - physics_variables.beta_poloidal_vol_avg = calculate_poloidal_beta( - physics_variables.b_plasma_total, - physics_variables.b_plasma_poloidal_average, - physics_variables.beta_total_vol_avg, + physics_variables.beta_poloidal_vol_avg = self.beta.calculate_poloidal_beta( + b_plasma_total=physics_variables.b_plasma_total, + b_plasma_poloidal_average=physics_variables.b_plasma_poloidal_average, + beta=physics_variables.beta_total_vol_avg, ) physics_variables.beta_thermal_vol_avg = ( @@ -8340,21 +8340,6 @@ def calculate_plasma_masses( ) -def calculate_poloidal_beta(b_plasma_total, b_plasma_poloidal_average, beta): - """Calculates total poloidal beta - - Author: James Morris (UKAEA) - - J.P. Freidberg, "Plasma physics and fusion energy", Cambridge University Press (2007) - Page 270 ISBN 0521851076 - - :param b_plasma_total: sum of the toroidal and poloidal fields (T) - :param b_plasma_poloidal_average: poloidal field (T) - :param beta: total plasma beta - """ - return beta * (b_plasma_total / b_plasma_poloidal_average) ** 2 - - def res_diff_time(rmajor, res_plasma, kappa95): """Calculates resistive diffusion time @@ -8841,6 +8826,27 @@ def calculate_beta_limit_from_norm( / (rminor * b_plasma_toroidal_on_axis) ) + @staticmethod + def calculate_poloidal_beta( + b_plasma_total: float, b_plasma_poloidal_average: float, beta: float + ) -> float: + """Calculates total poloidal beta (β_p) + + :type b_plasma_total: float + :param b_plasma_poloidal_average: The average poloidal magnetic field of the plasma (in Tesla). + :type b_plasma_poloidal_average: float + :param beta: The plasma beta, a dimensionless parameter representing the ratio of plasma pressure to magnetic pressure. + :type beta: float + :return: The calculated total poloidal beta. + :rtype: float + + :references: + - J.P. Freidberg, "Plasma physics and fusion energy", Cambridge University Press (2007) + Page 270 ISBN 0521851076 + + """ + return beta * (b_plasma_total / b_plasma_poloidal_average) ** 2 + @staticmethod def fast_alpha_beta( b_plasma_poloidal_average: float, diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index e34f0d715..3736b45c9 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -26,7 +26,6 @@ PlasmaBeta, calculate_current_coefficient_hastie, calculate_plasma_current_peng, - calculate_poloidal_beta, calculate_poloidal_field, calculate_volt_second_requirements, diamagnetic_fraction_hender, @@ -61,7 +60,7 @@ def physics(): def test_calculate_poloidal_beta(): """Test calculate_poloidal_beta()""" - beta_poloidal = calculate_poloidal_beta(5.347, 0.852, 0.0307) + beta_poloidal = PlasmaBeta.calculate_poloidal_beta(5.347, 0.852, 0.0307) assert beta_poloidal == pytest.approx(1.209, abs=0.001) From 347967185a55f37adccc3806b869285805839d2d Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 21 Jan 2026 09:35:08 +0000 Subject: [PATCH 11/23] Refactor thermal beta normalization calculation to use PlasmaBeta class method --- process/physics.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/process/physics.py b/process/physics.py index 545a9230c..facda606f 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1800,13 +1800,13 @@ def physics(self): # ======================================================= - physics_variables.beta_norm_thermal = ( - 1.0e8 - * physics_variables.beta_thermal_vol_avg - * physics_variables.rminor - * physics_variables.b_plasma_toroidal_on_axis - / physics_variables.plasma_current + physics_variables.beta_norm_thermal = self.beta.calculate_normalised_beta( + beta=physics_variables.beta_thermal_vol_avg, + rminor=physics_variables.rminor, + c_plasma=physics_variables.plasma_current, + b_field=physics_variables.b_plasma_toroidal_on_axis, ) + physics_variables.beta_norm_toroidal = ( physics_variables.beta_norm_total * ( From 434b971ae8c96ed3fa40f12b1355be3fd9e4c43b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 21 Jan 2026 09:54:24 +0000 Subject: [PATCH 12/23] Refactor beta norm max calculation to use PlasmaBeta class property for improved maintainability --- process/physics.py | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/process/physics.py b/process/physics.py index facda606f..9b0afe1ec 100644 --- a/process/physics.py +++ b/process/physics.py @@ -2651,19 +2651,9 @@ def physics(self): ) ) - # Map calculation methods to a dictionary - beta_norm_max_calculations = { - 0: physics_variables.beta_norm_max, - 1: physics_variables.beta_norm_max_wesson, - 2: physics_variables.beta_norm_max_original_scaling, - 3: physics_variables.beta_norm_max_menard, - 4: physics_variables.beta_norm_max_thloreus, - 5: physics_variables.beta_norm_max_stambaugh, - } - # Calculate beta_norm_max based on i_beta_norm_max - if int(physics_variables.i_beta_norm_max) in beta_norm_max_calculations: - physics_variables.beta_norm_max = beta_norm_max_calculations[ + if int(physics_variables.i_beta_norm_max) in self.beta.beta_norm_max_model_map: + physics_variables.beta_norm_max = self.beta.beta_norm_max_model_map[ int(physics_variables.i_beta_norm_max) ] else: @@ -8594,6 +8584,18 @@ def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE + @property + def beta_norm_max_model_map(self): + """Mapping of beta norm max model indices to their calculated values.""" + return { + 0: physics_variables.beta_norm_max, + 1: physics_variables.beta_norm_max_wesson, + 2: physics_variables.beta_norm_max_original_scaling, + 3: physics_variables.beta_norm_max_menard, + 4: physics_variables.beta_norm_max_thloreus, + 5: physics_variables.beta_norm_max_stambaugh, + } + @staticmethod def calculate_beta_norm_max_wesson(ind_plasma_internal_norm: float) -> float: """ From 8003b11e4a087ff7d08ebd76e9a9711f0c6fea22 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 21 Jan 2026 13:51:22 +0000 Subject: [PATCH 13/23] Refactor beta norm max calculation to use BetaNormMaxModel Enum method and improve error handling --- process/physics.py | 43 +++++++++++++++++++++++++++---------------- 1 file changed, 27 insertions(+), 16 deletions(-) diff --git a/process/physics.py b/process/physics.py index 9b0afe1ec..c308b50a1 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1,5 +1,6 @@ import logging import math +from enum import IntEnum import numba as nb import numpy as np @@ -2652,15 +2653,14 @@ def physics(self): ) # Calculate beta_norm_max based on i_beta_norm_max - if int(physics_variables.i_beta_norm_max) in self.beta.beta_norm_max_model_map: - physics_variables.beta_norm_max = self.beta.beta_norm_max_model_map[ - int(physics_variables.i_beta_norm_max) - ] - else: + try: + model = BetaNormMaxModel(int(physics_variables.i_beta_norm_max)) + physics_variables.beta_norm_max = self.beta.get_beta_norm_max_value(model) + except ValueError: raise ProcessValueError( "Illegal value of i_beta_norm_max", i_beta_norm_max=physics_variables.i_beta_norm_max, - ) + ) from None # calculate_beta_limit() returns the beta_vol_avg_max for beta physics_variables.beta_vol_avg_max = self.beta.calculate_beta_limit_from_norm( @@ -8577,6 +8577,17 @@ def reinke_tsep(b_plasma_toroidal_on_axis, flh, qstar, rmajor, eps, fgw, kappa, ) +class BetaNormMaxModel(IntEnum): + """Beta norm max (β_N_max) model types""" + + USER_INPUT = 0 + WESSON = 1 + ORIGINAL_SCALING = 2 + MENARD = 3 + THLOREUS = 4 + STAMBAUGH = 5 + + class PlasmaBeta: """Class to hold plasma beta calculations for plasma processing.""" @@ -8584,17 +8595,17 @@ def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE - @property - def beta_norm_max_model_map(self): - """Mapping of beta norm max model indices to their calculated values.""" - return { - 0: physics_variables.beta_norm_max, - 1: physics_variables.beta_norm_max_wesson, - 2: physics_variables.beta_norm_max_original_scaling, - 3: physics_variables.beta_norm_max_menard, - 4: physics_variables.beta_norm_max_thloreus, - 5: physics_variables.beta_norm_max_stambaugh, + def get_beta_norm_max_value(self, model: BetaNormMaxModel) -> float: + """Get the beta norm max value (β_N_max) for the specified model.""" + model_map = { + BetaNormMaxModel.USER_INPUT: physics_variables.beta_norm_max, + BetaNormMaxModel.WESSON: physics_variables.beta_norm_max_wesson, + BetaNormMaxModel.ORIGINAL_SCALING: physics_variables.beta_norm_max_original_scaling, + BetaNormMaxModel.MENARD: physics_variables.beta_norm_max_menard, + BetaNormMaxModel.THLOREUS: physics_variables.beta_norm_max_thloreus, + BetaNormMaxModel.STAMBAUGH: physics_variables.beta_norm_max_stambaugh, } + return model_map[model] @staticmethod def calculate_beta_norm_max_wesson(ind_plasma_internal_norm: float) -> float: From 09176294f4c0fb913ed9474db846261ccf581493 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 28 Jan 2026 09:26:36 +0000 Subject: [PATCH 14/23] Update plasma beta documentation to include function references for Wesson, Original Scaling, Menard, Tholerus, and Stambaugh relations --- .../physics-models/plasma_beta/plasma_beta.md | 23 ++++++++++++++----- 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/documentation/physics-models/plasma_beta/plasma_beta.md b/documentation/physics-models/plasma_beta/plasma_beta.md index b44fe0f15..286b260a5 100644 --- a/documentation/physics-models/plasma_beta/plasma_beta.md +++ b/documentation/physics-models/plasma_beta/plasma_beta.md @@ -170,7 +170,7 @@ beta_norm_max = 3.0 --------- -#### Wesson Relation +#### Wesson Relation | `calculate_beta_norm_max_wesson()` This can be activated by stating `i_beta_norm_max = 1` in the input file. @@ -188,7 +188,7 @@ This is only recommended for high aspect ratio tokamaks[^3]. --------- -#### Original Scaling Law +#### Original Scaling Law | `calculate_beta_norm_max_original()` This can be activated by stating `i_beta_norm_max = 2` in the input file. @@ -263,7 +263,7 @@ $$ --------- -#### Menard Beta Relation +#### Menard Beta Relation | `calculate_beta_norm_max_menard()` This can be activated by stating `i_beta_norm_max = 3` in the input file. @@ -345,7 +345,7 @@ Found as a reasonable fit to the computed no wall limit at $f_{\text{BS}} \appro --------- -#### Tholerus Relation +#### Tholerus Relation | `calculate_beta_norm_max_thloreus()` This can be activated by stating `i_beta_norm_max = 4` in the input file. @@ -362,7 +362,7 @@ where $F_p$ is the pressure peaking, $F_p = p_{\text{ax}} / \langle p \rangle$ a ------------- -#### Stambaugh Relation +#### Stambaugh Relation | `calculate_beta_norm_max_stambaugh()` This can be activated by stating `i_beta_norm_max = 5` in the input file. @@ -377,13 +377,24 @@ This fit was done for $A = 1.2 -7.0, \kappa = 1.5-6.0$ with $\delta = 0.5$ for n --------- +## Stored energy | `calculate_plasma_energy_from_beta()` + +As the $\beta$ metric is simply the ratio between the kinetic pressure of the plasma and the magnetic pressure, $\beta$ can be used to get the total kinetic energy of the plasma (assuming the total $\beta$ is used). + + +$$ +E_{\text{plasma}} \approx \frac{3}{2}\frac{\beta B^2}{2\mu_0}V_{\text{plasma}} +$$ + +--------------- + ## Key Constraints ### Beta consistency This constraint can be activated 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 $\mathtt{beta}$ is consistent in value with the rest of the physics parameters +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 $$ \texttt{beta_total_vol_avg} \equiv \frac{2\mu_0 \langle n_{\text{e}}T_{\text{e}}+n_{\text{i}}T_{\text{i}}\rangle}{B^2} + \beta_{\alpha} + \beta_{\text{beam}} From cab45228f7a74b66876a2a72d5827dbf15b083a6 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 28 Jan 2026 10:00:58 +0000 Subject: [PATCH 15/23] Refactor plasma current documentation for clarity and consistency in variable formatting --- .../plasma_current/plasma_current.md | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/documentation/physics-models/plasma_current/plasma_current.md b/documentation/physics-models/plasma_current/plasma_current.md index cce5fec99..fdd994fd9 100644 --- a/documentation/physics-models/plasma_current/plasma_current.md +++ b/documentation/physics-models/plasma_current/plasma_current.md @@ -92,7 +92,7 @@ Instead of $q_a$, $q_{95}$ is used as in plasma configurations with divertors th ## Plasma Current Calculation | `calculate_plasma_current()` -This function calculates the plasma current shaping factor ($f_q$), plasma current ($I_{\text{p}}$), qstar ($q^*$), normalized beta ($\beta_{\text{N}}$) and then poloidal field and the profile settings for $\mathtt{alphaj}$ ($\alpha_J$) and $\mathtt{ind_plasma_internal_norm}$ ($l_{\mathtt{i}}$). This is done in 5 separate steps which are shown in the following numbered sections. +This function calculates the plasma current shaping factor ($f_q$), plasma current ($I_{\text{p}}$), qstar ($q^*$) and then poloidal field and the profile settings for $\texttt{alphaj}$ ($\alpha_J$) and $\texttt{ind_plasma_internal_norm}$ ($l_{\text{i}}$). This is done in 5 separate steps which are shown in the following numbered sections. $$\begin{aligned} @@ -118,7 +118,7 @@ The formula for calculating `fq` is: $$f_q = \left(\frac{{1.22 - 0.68 \epsilon}}{{(1.0 - \epsilon^2)^2}}\right) \left(\frac{L_{\text{poloidal}}}{2\pi a}\right)^2$$ -Where $\epsilon$ is the inverse [aspect ratio](../plasma_geometry.md) ($\mathtt{eps}$) and $L_{\text{poloidal}}$ is the poloidal perimeter of the plasma. +Where $\epsilon$ is the inverse [aspect ratio](../plasma_geometry.md) ($\texttt{eps}$) and $L_{\text{poloidal}}$ is the poloidal perimeter of the plasma. ----------- @@ -523,23 +523,13 @@ $$ -------------- -### 3. Calculate the normalized beta -The total normalized beta is calculated as per: - -$$ -\beta_N = \beta\frac{1\times10^8 a B_{\text{T}}}{I_{\text{P}}} -$$ - - ------------------ - -### 4. Plasma Current Poloidal Field +### 3. Plasma Current Poloidal Field For calculating the poloidal magnetic field created due to the presence of the plasma current, [Ampere's law](https://en.wikipedia.org/wiki/Amp%C3%A8re%27s_circuital_law) can be used. In this case the poloidal field is simply returned as: $$ -B_{\text{p}} = \frac{\mu_0 I_{\text{p}}}{\mathtt{len_plasma_poloidal}} +B_{\text{p}} = \frac{\mu_0 I_{\text{p}}}{\texttt{len_plasma_poloidal}} $$ Where `len_plasma_poloidal` is the plasma poloidal perimeter calculated [here](../plasma_geometry.md#poloidal-perimeter). From c75741c3fb6473ed16ce9009b822fddadbe9ec54 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Wed, 4 Feb 2026 11:21:19 +0000 Subject: [PATCH 16/23] Refactor plasma beta calculations into a dedicated run method for improved organization and maintainability --- process/physics.py | 242 +++++++++++++++++++++++---------------------- 1 file changed, 125 insertions(+), 117 deletions(-) diff --git a/process/physics.py b/process/physics.py index c308b50a1..a1e582dde 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1743,107 +1743,7 @@ def physics(self): # Beta Components # ----------------------------------------------------- - physics_variables.beta_toroidal_vol_avg = ( - physics_variables.beta_total_vol_avg - * physics_variables.b_plasma_total**2 - / physics_variables.b_plasma_toroidal_on_axis**2 - ) - - # Calculate physics_variables.beta poloidal [-] - physics_variables.beta_poloidal_vol_avg = self.beta.calculate_poloidal_beta( - b_plasma_total=physics_variables.b_plasma_total, - b_plasma_poloidal_average=physics_variables.b_plasma_poloidal_average, - beta=physics_variables.beta_total_vol_avg, - ) - - physics_variables.beta_thermal_vol_avg = ( - physics_variables.beta_total_vol_avg - - physics_variables.beta_fast_alpha - - physics_variables.beta_beam - ) - - physics_variables.beta_poloidal_eps = ( - physics_variables.beta_poloidal_vol_avg * physics_variables.eps - ) - - physics_variables.beta_thermal_poloidal_vol_avg = ( - physics_variables.beta_thermal_vol_avg - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_poloidal_average - ) - ** 2 - ) - physics_variables.beta_thermal_toroidal_vol_avg = ( - physics_variables.beta_thermal_vol_avg - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_toroidal_on_axis - ) - ** 2 - ) - - # ======================================================= - - # Mirror the pressure profiles to match the doubled toroidal field profile - pres_profile_total = np.concatenate([ - physics_variables.pres_plasma_thermal_total_profile[::-1], - physics_variables.pres_plasma_thermal_total_profile, - ]) - - physics_variables.beta_thermal_toroidal_profile = np.array([ - self.calculate_plasma_beta( - pres_plasma=pres_profile_total[i], - b_field=physics_variables.b_plasma_toroidal_profile[i], - ) - for i in range(len(physics_variables.b_plasma_toroidal_profile)) - ]) - - # ======================================================= - - physics_variables.beta_norm_thermal = self.beta.calculate_normalised_beta( - beta=physics_variables.beta_thermal_vol_avg, - rminor=physics_variables.rminor, - c_plasma=physics_variables.plasma_current, - b_field=physics_variables.b_plasma_toroidal_on_axis, - ) - - physics_variables.beta_norm_toroidal = ( - physics_variables.beta_norm_total - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_toroidal_on_axis - ) - ** 2 - ) - physics_variables.beta_norm_poloidal = ( - physics_variables.beta_norm_total - * ( - physics_variables.b_plasma_total - / physics_variables.b_plasma_poloidal_average - ) - ** 2 - ) - - physics_variables.f_beta_alpha_beam_thermal = ( - physics_variables.beta_fast_alpha + physics_variables.beta_beam - ) / physics_variables.beta_thermal_vol_avg - - # Plasma thermal energy derived from the thermal beta - physics_variables.e_plasma_beta_thermal = ( - self.beta.calculate_plasma_energy_from_beta( - beta=physics_variables.beta_thermal_vol_avg, - b_field=physics_variables.b_plasma_total, - vol_plasma=physics_variables.vol_plasma, - ) - ) - - # Plasma thermal energy derived from the total beta - physics_variables.e_plasma_beta = self.beta.calculate_plasma_energy_from_beta( - beta=physics_variables.beta_total_vol_avg, - b_field=physics_variables.b_plasma_total, - vol_plasma=physics_variables.vol_plasma, - ) + self.beta.run() # ======================================================= @@ -3767,22 +3667,6 @@ def calculate_plasma_current( return b_plasma_poloidal_average, qstar, plasma_current - def calculate_plasma_beta(self, pres_plasma: float, b_field: float) -> float: - """ - Calculate the plasma beta for a given pressure and field. - - :param pres_plasma: Plasma pressure (in Pascals). - :type pres_plasma: float - :param b_field: Magnetic field strength (in Tesla). - :type b_field: float - :returns: The plasma beta (dimensionless). - :rtype: float - - Plasma beta is the ratio of plasma pressure to magnetic pressure. - """ - - return 2 * constants.RMU0 * pres_plasma / (b_field**2) - def outtim(self): po.oheadr(self.outfile, "Times") @@ -8607,6 +8491,130 @@ def get_beta_norm_max_value(self, model: BetaNormMaxModel) -> float: } return model_map[model] + def run(self) -> None: + """ + Calculate plasma beta values. + """ + + physics_variables.beta_toroidal_vol_avg = ( + physics_variables.beta_total_vol_avg + * physics_variables.b_plasma_total**2 + / physics_variables.b_plasma_toroidal_on_axis**2 + ) + + # Calculate physics_variables.beta poloidal [-] + physics_variables.beta_poloidal_vol_avg = self.calculate_poloidal_beta( + b_plasma_total=physics_variables.b_plasma_total, + b_plasma_poloidal_average=physics_variables.b_plasma_poloidal_average, + beta=physics_variables.beta_total_vol_avg, + ) + + physics_variables.beta_thermal_vol_avg = ( + physics_variables.beta_total_vol_avg + - physics_variables.beta_fast_alpha + - physics_variables.beta_beam + ) + + physics_variables.beta_poloidal_eps = ( + physics_variables.beta_poloidal_vol_avg * physics_variables.eps + ) + + physics_variables.beta_thermal_poloidal_vol_avg = ( + physics_variables.beta_thermal_vol_avg + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_poloidal_average + ) + ** 2 + ) + physics_variables.beta_thermal_toroidal_vol_avg = ( + physics_variables.beta_thermal_vol_avg + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_toroidal_on_axis + ) + ** 2 + ) + + # ======================================================= + + # Mirror the pressure profiles to match the doubled toroidal field profile + pres_profile_total = np.concatenate([ + physics_variables.pres_plasma_thermal_total_profile[::-1], + physics_variables.pres_plasma_thermal_total_profile, + ]) + + physics_variables.beta_thermal_toroidal_profile = np.array([ + self.calculate_plasma_beta( + pres_plasma=pres_profile_total[i], + b_field=physics_variables.b_plasma_toroidal_profile[i], + ) + for i in range(len(physics_variables.b_plasma_toroidal_profile)) + ]) + + # ======================================================= + + physics_variables.beta_norm_thermal = self.calculate_normalised_beta( + beta=physics_variables.beta_thermal_vol_avg, + rminor=physics_variables.rminor, + c_plasma=physics_variables.plasma_current, + b_field=physics_variables.b_plasma_toroidal_on_axis, + ) + + physics_variables.beta_norm_toroidal = ( + physics_variables.beta_norm_total + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_toroidal_on_axis + ) + ** 2 + ) + physics_variables.beta_norm_poloidal = ( + physics_variables.beta_norm_total + * ( + physics_variables.b_plasma_total + / physics_variables.b_plasma_poloidal_average + ) + ** 2 + ) + + physics_variables.f_beta_alpha_beam_thermal = ( + physics_variables.beta_fast_alpha + physics_variables.beta_beam + ) / physics_variables.beta_thermal_vol_avg + + # Plasma thermal energy derived from the thermal beta + physics_variables.e_plasma_beta_thermal = self.calculate_plasma_energy_from_beta( + beta=physics_variables.beta_thermal_vol_avg, + b_field=physics_variables.b_plasma_total, + vol_plasma=physics_variables.vol_plasma, + ) + + # Plasma thermal energy derived from the total beta + physics_variables.e_plasma_beta = self.calculate_plasma_energy_from_beta( + beta=physics_variables.beta_total_vol_avg, + b_field=physics_variables.b_plasma_total, + vol_plasma=physics_variables.vol_plasma, + ) + + @staticmethod + def calculate_plasma_beta( + pres_plasma: float | np.ndarray, b_field: float | np.ndarray + ) -> float | np.ndarray: + """ + Calculate the plasma beta (β) for a given pressure and field. + + :param pres_plasma: Plasma pressure (in Pascals). + :type pres_plasma: float | np.ndarray + :param b_field: Magnetic field strength (in Tesla). + :type b_field: float | np.ndarray + :returns: The plasma beta (dimensionless). + :rtype: float | np.ndarray + + Plasma beta is the ratio of plasma pressure to magnetic pressure. + """ + + return 2 * constants.RMU0 * pres_plasma / (b_field**2) + @staticmethod def calculate_beta_norm_max_wesson(ind_plasma_internal_norm: float) -> float: """ From 3ef1794455ecfc661816eec6d6f4e274cd38fd25 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 17:31:17 +0000 Subject: [PATCH 17/23] Add PlasmaInductance class and integrate into physics and stellarator models --- process/main.py | 5 +- process/physics.py | 270 +++++++++++++++++---------------- process/stellarator.py | 2 + tests/unit/test_physics.py | 5 +- tests/unit/test_stellarator.py | 4 +- 5 files changed, 151 insertions(+), 135 deletions(-) diff --git a/process/main.py b/process/main.py index 5d3460782..e9ca84254 100644 --- a/process/main.py +++ b/process/main.py @@ -94,7 +94,7 @@ ) from process.log import logging_model_handler, show_errors from process.pfcoil import PFCoil -from process.physics import DetailedPhysics, Physics, PlasmaBeta +from process.physics import DetailedPhysics, Physics, PlasmaBeta, PlasmaInductance from process.plasma_geometry import PlasmaGeom from process.plasma_profiles import PlasmaProfile from process.power import Power @@ -683,10 +683,12 @@ def __init__(self): electron_bernstein=ElectronBernstein(plasma_profile=self.plasma_profile), ) self.plasma_beta = PlasmaBeta() + self.plasma_inductance = PlasmaInductance() self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive, plasma_beta=self.plasma_beta, + plasma_inductance=self.plasma_inductance, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, @@ -704,6 +706,7 @@ def __init__(self): physics=self.physics, neoclassics=self.neoclassics, plasma_beta=self.plasma_beta, + plasma_inductance=self.plasma_inductance, ) self.dcll = DCLL(fw=self.fw) diff --git a/process/physics.py b/process/physics.py index a1e582dde..0d511474f 100644 --- a/process/physics.py +++ b/process/physics.py @@ -73,135 +73,6 @@ def rether( return conie * (temp_plasma_ion_vol_avg_kev - te) / (te**1.5) -@nb.jit(nopython=True, cache=True) -def calculate_volt_second_requirements( - csawth: float, - eps: float, - f_c_plasma_inductive: float, - ejima_coeff: float, - kappa: float, - rmajor: float, - res_plasma: float, - plasma_current: float, - t_plant_pulse_fusion_ramp: float, - t_plant_pulse_burn: float, - ind_plasma_internal_norm: float, -) -> tuple[float, float, float, float, float, float]: - """Calculate the volt-second requirements and related parameters for plasma physics. - - :param csawth: Coefficient for sawteeth effects - :type csawth: float - :param eps: Inverse aspect ratio - :type eps: float - :param f_c_plasma_inductive: Fraction of plasma current produced inductively - :type f_c_plasma_inductive: float - :param ejima_coeff: Ejima coefficient for resistive start-up V-s component - :type ejima_coeff: float - :param kappa: Plasma elongation - :type kappa: float - :param rmajor: Plasma major radius (m) - :type rmajor: float - :param res_plasma: Plasma resistance (ohm) - :type res_plasma: float - :param plasma_current: Plasma current (A) - :type plasma_current: float - :param t_plant_pulse_fusion_ramp: Heating time (s) - :type t_plant_pulse_fusion_ramp: float - :param t_plant_pulse_burn: Burn time (s) - :type t_plant_pulse_burn: float - :param ind_plasma_internal_norm: Plasma normalized internal inductance - :type ind_plasma_internal_norm: float - - :return: A tuple containing: - - vs_plasma_internal: Internal plasma volt-seconds (Wb) - - ind_plasma_internal: Plasma inductance (H) - - vs_plasma_burn_required: Volt-seconds needed during flat-top (heat+burn) (Wb) - - vs_plasma_ramp_required: Volt-seconds needed during ramp-up (Wb) - - ind_plasma_total,: Internal and external plasma inductance V-s (Wb) - - vs_res_ramp: Resistive losses in start-up volt-seconds (Wb) - - vs_plasma_total_required: Total volt-seconds needed (Wb) - :rtype: tuple[float, float, float, float, float, float] - - :notes: - - :references: - - S. Ejima, R. W. Callis, J. L. Luxon, R. D. Stambaugh, T. S. Taylor, and J. C. Wesley, - “Volt-second analysis and consumption in Doublet III plasmas,” - Nuclear Fusion, vol. 22, no. 10, pp. 1313-1319, Oct. 1982, doi: - https://doi.org/10.1088/0029-5515/22/10/006. - - - S. C. Jardin, C. E. Kessel, and N Pomphrey, - “Poloidal flux linkage requirements for the International Thermonuclear Experimental Reactor,” - Nuclear Fusion, vol. 34, no. 8, pp. 1145-1160, Aug. 1994, - doi: https://doi.org/10.1088/0029-5515/34/8/i07. - - - S. P. Hirshman and G. H. Neilson, “External inductance of an axisymmetric plasma,” - The Physics of Fluids, vol. 29, no. 3, pp. 790-793, Mar. 1986, - doi: https://doi.org/10.1063/1.865934. - """ - # Plasma internal inductance - - ind_plasma_internal = constants.RMU0 * rmajor * ind_plasma_internal_norm / 2.0 - - # Internal plasma flux (V-s) component - vs_plasma_internal = ind_plasma_internal * plasma_current - - # Start-up resistive component - # Uses ITER formula without the 10 V-s add-on - - vs_res_ramp = ejima_coeff * constants.RMU0 * plasma_current * rmajor - - # ====================================================================== - - # Hirshman and Neilson fit for external inductance - - aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( - 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps - ) - beps = 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) - - ind_plasma_external = ( - rmajor * constants.RMU0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) - ) - - # ====================================================================== - - ind_plasma_total = ind_plasma_external + ind_plasma_internal - - # Inductive V-s component - - vs_self_ind_ramp = ind_plasma_total * plasma_current - vs_plasma_ramp_required = vs_res_ramp + vs_self_ind_ramp - - # Plasma loop voltage during flat-top - # Include enhancement factor in flattop V-s requirement - # to account for MHD sawtooth effects. - - v_plasma_loop_burn = plasma_current * res_plasma * f_c_plasma_inductive - - v_burn_resistive = v_plasma_loop_burn * csawth - - # N.B. t_plant_pulse_burn on first iteration will not be correct - # if the pulsed reactor option is used, but the value - # will be correct on subsequent calls. - - vs_plasma_burn_required = v_burn_resistive * ( - t_plant_pulse_fusion_ramp + t_plant_pulse_burn - ) - vs_plasma_total_required = vs_plasma_ramp_required + vs_plasma_burn_required - - return ( - vs_plasma_internal, - ind_plasma_total, - vs_plasma_burn_required, - vs_plasma_ramp_required, - vs_self_ind_ramp, - vs_res_ramp, - vs_plasma_total_required, - v_plasma_loop_burn, - ) - - # ----------------------------------------------------- # Plasma Current & Poloidal Field Calculations # ----------------------------------------------------- @@ -1516,12 +1387,13 @@ def _trapped_particle_fraction_sauter( class Physics: - def __init__(self, plasma_profile, current_drive, plasma_beta): + def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance): self.outfile = constants.NOUT self.mfile = constants.MFILE self.plasma_profile = plasma_profile self.current_drive = current_drive self.beta = plasma_beta + self.inductance = plasma_inductance def physics(self): """ @@ -2460,7 +2332,7 @@ def physics(self): physics_variables.vs_plasma_res_ramp, physics_variables.vs_plasma_total_required, physics_variables.v_plasma_loop_burn, - ) = calculate_volt_second_requirements( + ) = self.inductance.calculate_volt_second_requirements( physics_variables.csawth, physics_variables.eps, physics_variables.f_c_plasma_inductive, @@ -9237,6 +9109,142 @@ def output_beta_information(self): po.oblnkl(self.outfile) +class PlasmaInductance: + """Class to hold plasma inductance calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + @staticmethod + def calculate_volt_second_requirements( + csawth: float, + eps: float, + f_c_plasma_inductive: float, + ejima_coeff: float, + kappa: float, + rmajor: float, + res_plasma: float, + plasma_current: float, + t_plant_pulse_fusion_ramp: float, + t_plant_pulse_burn: float, + ind_plasma_internal_norm: float, + ) -> tuple[float, float, float, float, float, float]: + """Calculate the volt-second requirements and related parameters for plasma physics. + + :param csawth: Coefficient for sawteeth effects + :type csawth: float + :param eps: Inverse aspect ratio + :type eps: float + :param f_c_plasma_inductive: Fraction of plasma current produced inductively + :type f_c_plasma_inductive: float + :param ejima_coeff: Ejima coefficient for resistive start-up V-s component + :type ejima_coeff: float + :param kappa: Plasma elongation + :type kappa: float + :param rmajor: Plasma major radius (m) + :type rmajor: float + :param res_plasma: Plasma resistance (ohm) + :type res_plasma: float + :param plasma_current: Plasma current (A) + :type plasma_current: float + :param t_plant_pulse_fusion_ramp: Heating time (s) + :type t_plant_pulse_fusion_ramp: float + :param t_plant_pulse_burn: Burn time (s) + :type t_plant_pulse_burn: float + :param ind_plasma_internal_norm: Plasma normalized internal inductance + :type ind_plasma_internal_norm: float + + :return: A tuple containing: + - vs_plasma_internal: Internal plasma volt-seconds (Wb) + - ind_plasma_internal: Plasma inductance (H) + - vs_plasma_burn_required: Volt-seconds needed during flat-top (heat+burn) (Wb) + - vs_plasma_ramp_required: Volt-seconds needed during ramp-up (Wb) + - ind_plasma_total,: Internal and external plasma inductance V-s (Wb) + - vs_res_ramp: Resistive losses in start-up volt-seconds (Wb) + - vs_plasma_total_required: Total volt-seconds needed (Wb) + :rtype: tuple[float, float, float, float, float, float] + + :notes: + + :references: + - S. Ejima, R. W. Callis, J. L. Luxon, R. D. Stambaugh, T. S. Taylor, and J. C. Wesley, + “Volt-second analysis and consumption in Doublet III plasmas,” + Nuclear Fusion, vol. 22, no. 10, pp. 1313-1319, Oct. 1982, doi: + https://doi.org/10.1088/0029-5515/22/10/006. + + - S. C. Jardin, C. E. Kessel, and N Pomphrey, + “Poloidal flux linkage requirements for the International Thermonuclear Experimental Reactor,” + Nuclear Fusion, vol. 34, no. 8, pp. 1145-1160, Aug. 1994, + doi: https://doi.org/10.1088/0029-5515/34/8/i07. + + - S. P. Hirshman and G. H. Neilson, “External inductance of an axisymmetric plasma,” + The Physics of Fluids, vol. 29, no. 3, pp. 790-793, Mar. 1986, + doi: https://doi.org/10.1063/1.865934. + """ + # Plasma internal inductance + + ind_plasma_internal = constants.RMU0 * rmajor * ind_plasma_internal_norm / 2.0 + + # Internal plasma flux (V-s) component + vs_plasma_internal = ind_plasma_internal * plasma_current + + # Start-up resistive component + # Uses ITER formula without the 10 V-s add-on + + vs_res_ramp = ejima_coeff * constants.RMU0 * plasma_current * rmajor + + # ====================================================================== + + # Hirshman and Neilson fit for external inductance + + aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( + 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps + ) + beps = 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) + + ind_plasma_external = ( + rmajor * constants.RMU0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) + ) + + # ====================================================================== + + ind_plasma_total = ind_plasma_external + ind_plasma_internal + + # Inductive V-s component + + vs_self_ind_ramp = ind_plasma_total * plasma_current + vs_plasma_ramp_required = vs_res_ramp + vs_self_ind_ramp + + # Plasma loop voltage during flat-top + # Include enhancement factor in flattop V-s requirement + # to account for MHD sawtooth effects. + + v_plasma_loop_burn = plasma_current * res_plasma * f_c_plasma_inductive + + v_burn_resistive = v_plasma_loop_burn * csawth + + # N.B. t_plant_pulse_burn on first iteration will not be correct + # if the pulsed reactor option is used, but the value + # will be correct on subsequent calls. + + vs_plasma_burn_required = v_burn_resistive * ( + t_plant_pulse_fusion_ramp + t_plant_pulse_burn + ) + vs_plasma_total_required = vs_plasma_ramp_required + vs_plasma_burn_required + + return ( + vs_plasma_internal, + ind_plasma_total, + vs_plasma_burn_required, + vs_plasma_ramp_required, + vs_self_ind_ramp, + vs_res_ramp, + vs_plasma_total_required, + v_plasma_loop_burn, + ) + + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" diff --git a/process/stellarator.py b/process/stellarator.py index 18533d959..e24ff0d58 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -67,6 +67,7 @@ def __init__( physics, neoclassics, plasma_beta, + plasma_inductance, ) -> None: """Initialises the Stellarator model's variables @@ -104,6 +105,7 @@ def __init__( self.physics = physics self.neoclassics = neoclassics self.beta = plasma_beta + self.inductance = plasma_inductance def run(self, output: bool): """Routine to call the physics and engineering modules diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 3736b45c9..3efcb2c05 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -24,10 +24,10 @@ DetailedPhysics, Physics, PlasmaBeta, + PlasmaInductance, calculate_current_coefficient_hastie, calculate_plasma_current_peng, calculate_poloidal_field, - calculate_volt_second_requirements, diamagnetic_fraction_hender, diamagnetic_fraction_scene, ps_fraction_scene, @@ -55,6 +55,7 @@ def physics(): lower_hybrid=LowerHybrid(plasma_profile=PlasmaProfile()), ), PlasmaBeta(), + PlasmaInductance(), ) @@ -2099,7 +2100,7 @@ def test_vscalc(voltsecondreqparam): vs_plasma_res_ramp, vs_plasma_total_required, v_plasma_loop_burn, - ) = calculate_volt_second_requirements( + ) = PlasmaInductance.calculate_volt_second_requirements( csawth=voltsecondreqparam.csawth, eps=voltsecondreqparam.eps, f_c_plasma_inductive=voltsecondreqparam.f_c_plasma_inductive, diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 4d3320414..b72371607 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -28,7 +28,7 @@ ) from process.fw import Fw from process.hcpb import CCFE_HCPB -from process.physics import Physics, PlasmaBeta +from process.physics import Physics, PlasmaBeta, PlasmaInductance from process.plasma_profiles import PlasmaProfile from process.power import Power from process.stellarator import Neoclassics, Stellarator @@ -69,9 +69,11 @@ def stellarator(): ElectronBernstein(plasma_profile=PlasmaProfile()), ), PlasmaBeta(), + PlasmaInductance(), ), Neoclassics(), plasma_beta=PlasmaBeta(), + plasma_inductance=PlasmaInductance(), ) From 6bfea611869313604a7f16db9a072de34d0d04af Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 18:01:11 +0000 Subject: [PATCH 18/23] Refactor internal inductance calculations to use PlasmaInductance class methods and remove redundant static methods --- process/physics.py | 180 +++++++++++++++++++++++---------------------- 1 file changed, 92 insertions(+), 88 deletions(-) diff --git a/process/physics.py b/process/physics.py index 0d511474f..e6c5666dc 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1506,17 +1506,21 @@ def physics(self): # ----------------------------------------------------- physics_variables.ind_plasma_internal_norm_wesson = ( - self.calculate_internal_inductance_wesson(alphaj=physics_variables.alphaj) + self.inductance.calculate_internal_inductance_wesson( + alphaj=physics_variables.alphaj + ) ) # Spherical Tokamak relation for internal inductance # Menard et al. (2016), Nuclear Fusion, 56, 106023 physics_variables.ind_plasma_internal_norm_menard = ( - self.calculate_internal_inductance_menard(kappa=physics_variables.kappa) + self.inductance.calculate_internal_inductance_menard( + kappa=physics_variables.kappa + ) ) physics_variables.ind_plasma_internal_norm_iter_3 = ( - self.calculate_normalised_internal_inductance_iter_3( + self.inductance.calculate_normalised_internal_inductance_iter_3( b_plasma_poloidal_vol_avg=physics_variables.b_plasma_poloidal_average, c_plasma=physics_variables.plasma_current, vol_plasma=physics_variables.vol_plasma, @@ -2650,50 +2654,6 @@ def calculate_current_profile_index_wesson(qstar: float, q0: float) -> float: """ return qstar / q0 - 1.0 - @staticmethod - def calculate_internal_inductance_wesson(alphaj: float) -> float: - """ - Calculate the Wesson plasma normalized internal inductance. - - :param alphaj: Current profile index. - :type alphaj: float - - :return: The Wesson plasma normalised internal inductance. - :rtype: float - - :Notes: - - It is recommended to use this method with the other Wesson relations for normalised beta and - current profile index. - - This relation is only true for the cyclindrical plasma approximation with parabolic profiles. - - :References: - - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, - International Series of Monographs on Physics, Volume 149. - """ - return np.log(1.65 + 0.89 * alphaj) - - @staticmethod - def calculate_internal_inductance_menard(kappa: float) -> float: - """ - Calculate the Menard plasma normalized internal inductance. - - :param kappa: Plasma separatrix elongation. - :type kappa: float - - :return: The Menard plasma normalised internal inductance. - :rtype: float - - :Notes: - - This relation is based off of data from NSTX for l_i in the range of 0.4-0.85 - - This model is only recommneded to be used for ST's with kappa > 2.5 - - :References: - - J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” - Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, - doi: https://doi.org/10.1088/0029-5515/56/10/106023. - """ - return 3.4 - kappa - @staticmethod def calculate_density_limit( b_plasma_toroidal_on_axis: float, @@ -3254,47 +3214,6 @@ def phyaux( f_alpha_energy_confinement, ) - @staticmethod - def calculate_normalised_internal_inductance_iter_3( - b_plasma_poloidal_vol_avg: float, - c_plasma: float, - vol_plasma: float, - rmajor: float, - ) -> float: - """ - Calculate the normalised internal inductance using ITER-3 scaling li(3). - - :param b_plasma_poloidal_vol_avg: Volume-averaged poloidal magnetic field (T). - :type b_plasma_poloidal_vol_avg: float - :param c_plasma: Plasma current (A). - :type c_plasma: float - :param vol_plasma: Plasma volume (m^3). - :type vol_plasma: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - - :returns: The li(3) normalised internal inductance. - :rtype: float - - :references: - - T. C. Luce, D. A. Humphreys, G. L. Jackson, and W. M. Solomon, - “Inductive flux usage and its optimization in tokamak operation,” - Nuclear Fusion, vol. 54, no. 9, p. 093005, Jul. 2014, - doi: https://doi.org/10.1088/0029-5515/54/9/093005. - - - G. L. Jackson et al., “ITER startup studies in the DIII-D tokamak,” - Nuclear Fusion, vol. 48, no. 12, p. 125002, Nov. 2008, - doi: https://doi.org/10.1088/0029-5515/48/12/125002. - ‌ - """ - - return ( - 2 - * vol_plasma - * b_plasma_poloidal_vol_avg**2 - / (constants.RMU0**2 * c_plasma**2 * rmajor) - ) - @staticmethod def plasma_ohmic_heating( f_c_plasma_inductive: float, @@ -9244,6 +9163,91 @@ def calculate_volt_second_requirements( v_plasma_loop_burn, ) + @staticmethod + def calculate_normalised_internal_inductance_iter_3( + b_plasma_poloidal_vol_avg: float, + c_plasma: float, + vol_plasma: float, + rmajor: float, + ) -> float: + """ + Calculate the normalised internal inductance using ITER-3 scaling li(3). + + :param b_plasma_poloidal_vol_avg: Volume-averaged poloidal magnetic field (T). + :type b_plasma_poloidal_vol_avg: float + :param c_plasma: Plasma current (A). + :type c_plasma: float + :param vol_plasma: Plasma volume (m^3). + :type vol_plasma: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + + :returns: The li(3) normalised internal inductance. + :rtype: float + + :references: + - T. C. Luce, D. A. Humphreys, G. L. Jackson, and W. M. Solomon, + “Inductive flux usage and its optimization in tokamak operation,” + Nuclear Fusion, vol. 54, no. 9, p. 093005, Jul. 2014, + doi: https://doi.org/10.1088/0029-5515/54/9/093005. + + - G. L. Jackson et al., “ITER startup studies in the DIII-D tokamak,” + Nuclear Fusion, vol. 48, no. 12, p. 125002, Nov. 2008, + doi: https://doi.org/10.1088/0029-5515/48/12/125002. + ‌ + """ + + return ( + 2 + * vol_plasma + * b_plasma_poloidal_vol_avg**2 + / (constants.RMU0**2 * c_plasma**2 * rmajor) + ) + + @staticmethod + def calculate_internal_inductance_menard(kappa: float) -> float: + """ + Calculate the Menard plasma normalized internal inductance. + + :param kappa: Plasma separatrix elongation. + :type kappa: float + + :return: The Menard plasma normalised internal inductance. + :rtype: float + + :Notes: + - This relation is based off of data from NSTX for l_i in the range of 0.4-0.85 + - This model is only recommneded to be used for ST's with kappa > 2.5 + + :References: + - J. E. Menard et al., “Fusion nuclear science facilities and pilot plants based on the spherical tokamak,” + Nuclear Fusion, vol. 56, no. 10, p. 106023, Aug. 2016, + doi: https://doi.org/10.1088/0029-5515/56/10/106023. + """ + return 3.4 - kappa + + @staticmethod + def calculate_internal_inductance_wesson(alphaj: float) -> float: + """ + Calculate the Wesson plasma normalized internal inductance. + + :param alphaj: Current profile index. + :type alphaj: float + + :return: The Wesson plasma normalised internal inductance. + :rtype: float + + :Notes: + - It is recommended to use this method with the other Wesson relations for normalised beta and + current profile index. + - This relation is only true for the cyclindrical plasma approximation with parabolic profiles. + + :References: + - Wesson, J. (2011) Tokamaks. 4th Edition, 2011 Oxford Science Publications, + International Series of Monographs on Physics, Volume 149. + """ + return np.log(1.65 + 0.89 * alphaj) + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" From 01a860181d8722a078e149b61eaabf874e0bf0a6 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 18:09:32 +0000 Subject: [PATCH 19/23] Refactor internal inductance calculation to use IndInternalNormModel and PlasmaInductance class for improved clarity and maintainability --- process/physics.py | 42 +++++++++++++++++++++++--------------- tests/unit/test_physics.py | 6 +++--- 2 files changed, 28 insertions(+), 20 deletions(-) diff --git a/process/physics.py b/process/physics.py index e6c5666dc..92e3f2537 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1528,28 +1528,19 @@ def physics(self): ) ) - # Map calculation methods to a dictionary - ind_plasma_internal_norm_calculations = { - 0: physics_variables.ind_plasma_internal_norm, - 1: physics_variables.ind_plasma_internal_norm_wesson, - 2: physics_variables.ind_plasma_internal_norm_menard, - } - - # Calculate ind_plasma_internal_normx based on i_ind_plasma_internal_norm - if ( - int(physics_variables.i_ind_plasma_internal_norm) - in ind_plasma_internal_norm_calculations - ): + # Calculate ind_plasma_internal_norm based on i_ind_plasma_internal_norm + try: + model = IndInternalNormModel( + int(physics_variables.i_ind_plasma_internal_norm) + ) physics_variables.ind_plasma_internal_norm = ( - ind_plasma_internal_norm_calculations[ - int(physics_variables.i_ind_plasma_internal_norm) - ] + self.inductance.get_ind_internal_norm_value(model) ) - else: + except ValueError: raise ProcessValueError( "Illegal value of i_ind_plasma_internal_norm", i_ind_plasma_internal_norm=physics_variables.i_ind_plasma_internal_norm, - ) + ) from None # =================================================== @@ -9028,6 +9019,14 @@ def output_beta_information(self): po.oblnkl(self.outfile) +class IndInternalNormModel(IntEnum): + """Normalised internal inductance (l_i) model types""" + + USER_INPUT = 0 + WESSON = 1 + MENARD = 2 + + class PlasmaInductance: """Class to hold plasma inductance calculations for plasma processing.""" @@ -9035,6 +9034,15 @@ def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE + def get_ind_internal_norm_value(self, model: IndInternalNormModel) -> float: + """Get the normalised internal inductance (l_i) for the specified model.""" + model_map = { + IndInternalNormModel.USER_INPUT: physics_variables.ind_plasma_internal_norm, + IndInternalNormModel.WESSON: physics_variables.ind_plasma_internal_norm_wesson, + IndInternalNormModel.MENARD: physics_variables.ind_plasma_internal_norm_menard, + } + return model_map[model] + @staticmethod def calculate_volt_second_requirements( csawth: float, diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 3efcb2c05..9992a24a7 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -3400,14 +3400,14 @@ def test_calculate_current_profile_index_wesson(): def test_calculate_internal_inductance_wesson(): """Test calculate_internal_inductance_wesson().""" alphaj = 0.8 - result = Physics.calculate_internal_inductance_wesson(alphaj) + result = PlasmaInductance.calculate_internal_inductance_wesson(alphaj) assert result == pytest.approx(0.8595087177751706, abs=0.0001) def test_calculate_internal_inductance_menard(): """Test calculate_internal_inductance_menard().""" kappa = 2.8 - result = Physics.calculate_internal_inductance_menard(kappa) + result = PlasmaInductance.calculate_internal_inductance_menard(kappa) assert result == pytest.approx(0.6, abs=0.001) @@ -3456,7 +3456,7 @@ def test_calculate_beta_norm_max_stambaugh(): def test_calculate_internal_inductance_iter_3(): """Test calculate_normalised_internal_inductance_iter_3.""" - result = Physics.calculate_normalised_internal_inductance_iter_3( + result = PlasmaInductance.calculate_normalised_internal_inductance_iter_3( b_plasma_poloidal_vol_avg=1.0, c_plasma=1.5e7, vol_plasma=1000.0, rmajor=6.2 ) assert result == pytest.approx(0.9078959099585583, abs=0.00001) From 3a65b92ce7701c2323610824d80d55db89d640b3 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 20:38:45 +0000 Subject: [PATCH 20/23] Refactor output of volt-second information into PlasmaInductance class for improved organization and maintainability --- process/physics.py | 217 +++++++++++++++++++++++---------------------- 1 file changed, 111 insertions(+), 106 deletions(-) diff --git a/process/physics.py b/process/physics.py index 92e3f2537..a23b0208e 100644 --- a/process/physics.py +++ b/process/physics.py @@ -5611,112 +5611,7 @@ def outplas(self): po.ostars(self.outfile, 110) po.oblnkl(self.outfile) - po.osubhd(self.outfile, "Plasma Volt-second Requirements :") - po.ovarre( - self.outfile, - "Total plasma volt-seconds required for pulse (Wb)", - "(vs_plasma_total_required)", - physics_variables.vs_plasma_total_required, - "OP ", - ) - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Total plasma inductive flux consumption for plasma current ramp-up (Wb)", - "(vs_plasma_ind_ramp)", - physics_variables.vs_plasma_ind_ramp, - "OP ", - ) - po.ovarre( - self.outfile, - "Plasma resistive flux consumption for plasma current ramp-up (Wb)", - "(vs_plasma_res_ramp)", - physics_variables.vs_plasma_res_ramp, - "OP ", - ) - po.ovarre( - self.outfile, - "Total flux consumption for plasma current ramp-up (Wb)", - "(vs_plasma_ramp_required)", - physics_variables.vs_plasma_ramp_required, - "OP ", - ) - po.ovarrf( - self.outfile, - "Ejima coefficient", - "(ejima_coeff)", - physics_variables.ejima_coeff, - ) - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Internal plasma V-s", - "(vs_plasma_internal)", - physics_variables.vs_plasma_internal, - ) - - po.ovarre( - self.outfile, - "Plasma volt-seconds needed for flat-top (heat + burn times) (Wb)", - "(vs_plasma_burn_required)", - physics_variables.vs_plasma_burn_required, - "OP ", - ) - - po.ovarre( - self.outfile, - "Plasma loop voltage during burn (V)", - "(v_plasma_loop_burn)", - physics_variables.v_plasma_loop_burn, - "OP ", - ) - po.ovarrf( - self.outfile, - "Coefficient for sawtooth effects on burn V-s requirement", - "(csawth)", - physics_variables.csawth, - ) - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Plasma resistance (ohm)", - "(res_plasma)", - physics_variables.res_plasma, - "OP ", - ) - - po.ovarre( - self.outfile, - "Plasma resistive diffusion time (s)", - "(t_plasma_res_diffusion)", - physics_variables.t_plasma_res_diffusion, - "OP ", - ) - po.ovarre( - self.outfile, - "Plasma inductance (H)", - "(ind_plasma)", - physics_variables.ind_plasma, - "OP ", - ) - po.ovarre( - self.outfile, - "Plasma magnetic energy stored (J)", - "(e_plasma_magnetic_stored)", - physics_variables.e_plasma_magnetic_stored, - "OP ", - ) - po.ovarrf( - self.outfile, - "Plasma normalised internal inductance", - "(ind_plasma_internal_norm)", - physics_variables.ind_plasma_internal_norm, - "OP ", - ) - - po.oblnkl(self.outfile) - po.ostars(self.outfile, 110) - po.oblnkl(self.outfile) + self.inductance.output_volt_second_information() po.ovarrf( self.outfile, @@ -9256,6 +9151,116 @@ def calculate_internal_inductance_wesson(alphaj: float) -> float: """ return np.log(1.65 + 0.89 * alphaj) + def output_volt_second_information(self): + """Output volt-second information to file.""" + + po.osubhd(self.outfile, "Plasma Volt-second Requirements :") + po.ovarre( + self.outfile, + "Total plasma volt-seconds required for pulse (Wb)", + "(vs_plasma_total_required)", + physics_variables.vs_plasma_total_required, + "OP ", + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Total plasma inductive flux consumption for plasma current ramp-up (Wb)", + "(vs_plasma_ind_ramp)", + physics_variables.vs_plasma_ind_ramp, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma resistive flux consumption for plasma current ramp-up (Wb)", + "(vs_plasma_res_ramp)", + physics_variables.vs_plasma_res_ramp, + "OP ", + ) + po.ovarre( + self.outfile, + "Total flux consumption for plasma current ramp-up (Wb)", + "(vs_plasma_ramp_required)", + physics_variables.vs_plasma_ramp_required, + "OP ", + ) + po.ovarrf( + self.outfile, + "Ejima coefficient", + "(ejima_coeff)", + physics_variables.ejima_coeff, + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Internal plasma V-s", + "(vs_plasma_internal)", + physics_variables.vs_plasma_internal, + ) + + po.ovarre( + self.outfile, + "Plasma volt-seconds needed for flat-top (heat + burn times) (Wb)", + "(vs_plasma_burn_required)", + physics_variables.vs_plasma_burn_required, + "OP ", + ) + + po.ovarre( + self.outfile, + "Plasma loop voltage during burn (V)", + "(v_plasma_loop_burn)", + physics_variables.v_plasma_loop_burn, + "OP ", + ) + po.ovarrf( + self.outfile, + "Coefficient for sawtooth effects on burn V-s requirement", + "(csawth)", + physics_variables.csawth, + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Plasma resistance (ohm)", + "(res_plasma)", + physics_variables.res_plasma, + "OP ", + ) + + po.ovarre( + self.outfile, + "Plasma resistive diffusion time (s)", + "(t_plasma_res_diffusion)", + physics_variables.t_plasma_res_diffusion, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma inductance (H)", + "(ind_plasma)", + physics_variables.ind_plasma, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma magnetic energy stored (J)", + "(e_plasma_magnetic_stored)", + physics_variables.e_plasma_magnetic_stored, + "OP ", + ) + po.ovarrf( + self.outfile, + "Plasma normalised internal inductance", + "(ind_plasma_internal_norm)", + physics_variables.ind_plasma_internal_norm, + "OP ", + ) + + po.oblnkl(self.outfile) + po.ostars(self.outfile, 110) + po.oblnkl(self.outfile) + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" From ad631d10c4fc751a9e694998e6ee82476204c7f0 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 20:44:14 +0000 Subject: [PATCH 21/23] Update documentation for PlasmaInductance: add method references to Wesson and Menard relations --- .../physics-models/plasma_current/plasma_inductance.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/documentation/physics-models/plasma_current/plasma_inductance.md b/documentation/physics-models/plasma_current/plasma_inductance.md index 8644e6de0..0d24f321a 100644 --- a/documentation/physics-models/plasma_current/plasma_inductance.md +++ b/documentation/physics-models/plasma_current/plasma_inductance.md @@ -1,4 +1,4 @@ -# Plasma Inductance +# Plasma Inductance | `PlasmaInductance` ## Setting the normalised internal inductance @@ -22,7 +22,7 @@ ind_plasma_internal_norm = 1.0 ---------- -### Wesson relation +### Wesson relation | `calculate_internal_inductance_wesson()` This can be activated by stating `i_ind_plasma_internal_normx = 1` in the input file. @@ -42,7 +42,7 @@ This is only recommended for high aspect ratio tokamaks[^2]. --------- -### Menard Inductance Relation +### Menard Inductance Relation | `calculate_internal_inductance_menard()` This can be activated by stating `ind_plasma_internal_norm = 2` in the input file. From 9bcc89a79044d3a6c03536965d0cca7f2c4413b5 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 21:31:56 +0000 Subject: [PATCH 22/23] Add PlasmaCurrent class and integrate into Physics and Stellarator models --- process/main.py | 11 ++++++++++- process/physics.py | 18 +++++++++++++++++- process/stellarator.py | 2 ++ 3 files changed, 29 insertions(+), 2 deletions(-) diff --git a/process/main.py b/process/main.py index e9ca84254..f44dfe071 100644 --- a/process/main.py +++ b/process/main.py @@ -94,7 +94,13 @@ ) from process.log import logging_model_handler, show_errors from process.pfcoil import PFCoil -from process.physics import DetailedPhysics, Physics, PlasmaBeta, PlasmaInductance +from process.physics import ( + DetailedPhysics, + Physics, + PlasmaBeta, + PlasmaCurrent, + PlasmaInductance, +) from process.plasma_geometry import PlasmaGeom from process.plasma_profiles import PlasmaProfile from process.power import Power @@ -684,11 +690,13 @@ def __init__(self): ) self.plasma_beta = PlasmaBeta() self.plasma_inductance = PlasmaInductance() + self.plasma_current = PlasmaCurrent() self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive, plasma_beta=self.plasma_beta, plasma_inductance=self.plasma_inductance, + plasma_current=self.plasma_current, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, @@ -707,6 +715,7 @@ def __init__(self): neoclassics=self.neoclassics, plasma_beta=self.plasma_beta, plasma_inductance=self.plasma_inductance, + plasma_current=self.plasma_current, ) self.dcll = DCLL(fw=self.fw) diff --git a/process/physics.py b/process/physics.py index a23b0208e..1cec77f7b 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1387,13 +1387,21 @@ def _trapped_particle_fraction_sauter( class Physics: - def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance): + def __init__( + self, + plasma_profile, + current_drive, + plasma_beta, + plasma_inductance, + plasma_current, + ): self.outfile = constants.NOUT self.mfile = constants.MFILE self.plasma_profile = plasma_profile self.current_drive = current_drive self.beta = plasma_beta self.inductance = plasma_inductance + self.current = plasma_current def physics(self): """ @@ -9262,6 +9270,14 @@ def output_volt_second_information(self): po.oblnkl(self.outfile) +class PlasmaCurrent: + """Class to hold plasma current calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" diff --git a/process/stellarator.py b/process/stellarator.py index e24ff0d58..4f5a76054 100644 --- a/process/stellarator.py +++ b/process/stellarator.py @@ -68,6 +68,7 @@ def __init__( neoclassics, plasma_beta, plasma_inductance, + plasma_current, ) -> None: """Initialises the Stellarator model's variables @@ -106,6 +107,7 @@ def __init__( self.neoclassics = neoclassics self.beta = plasma_beta self.inductance = plasma_inductance + self.current = plasma_current def run(self, output: bool): """Routine to call the physics and engineering modules From 935e8f0975d57f859bec9fb094ea787e87845ca6 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 22:08:39 +0000 Subject: [PATCH 23/23] Refactor plasma current calculations to use PlasmaCurrent class - Updated test_physics.py to instantiate PlasmaCurrent for calculating plasma current and related methods. - Replaced direct calls to calculate_plasma_current_peng and calculate_current_coefficient_hastie with PlasmaCurrent class methods. - Modified test_stellarator.py to include PlasmaCurrent in the stellarator setup. --- process/physics.py | 1038 ++++++++++++++++---------------- tests/unit/test_physics.py | 42 +- tests/unit/test_stellarator.py | 4 +- 3 files changed, 547 insertions(+), 537 deletions(-) diff --git a/process/physics.py b/process/physics.py index 1cec77f7b..885192321 100644 --- a/process/physics.py +++ b/process/physics.py @@ -188,293 +188,6 @@ def calculate_poloidal_field( return b_plasma_toroidal_on_axis * (ff1 + ff2) / (2.0 * np.pi * qbar) -def calculate_current_coefficient_peng( - eps: float, len_plasma_poloidal: float, rminor: float -) -> float: - """ - Calculate the plasma current scaling coefficient for the Peng scaling from the STAR code. - - :param eps: Plasma inverse aspect ratio. - :type eps: float - :param len_plasma_poloidal: Plasma poloidal perimeter length [m]. - :type len_plasma_poloidal: float - :param rminor: Plasma minor radius [m]. - :type rminor: float - - :return: The plasma current scaling coefficient. - :rtype: float - - :references: None - """ - - return ( - (1.22 - 0.68 * eps) - / ((1.0 - eps * eps) ** 2) - * (len_plasma_poloidal / (2.0 * np.pi * rminor)) ** 2 - ) - - -def calculate_plasma_current_peng( - q95: float, - aspect: float, - eps: float, - rminor: float, - b_plasma_toroidal_on_axis: float, - kappa: float, - delta: float, -) -> float: - """ - Function to calculate plasma current (Peng scaling from the STAR code) - - Parameters: - - q95: float, 95% flux surface safety factor - - aspect: float, plasma aspect ratio - - eps: float, inverse aspect ratio - - rminor: float, plasma minor radius (m) - - b_plasma_toroidal_on_axis: float, toroidal field on axis (T) - - kappa: float, plasma elongation - - delta: float, plasma triangularity - - Returns: - - float, plasma current in MA - - This function calculates the plasma current in MA, - using a scaling from Peng, Galambos and Shipe (1992). - It is primarily used for Tight Aspect Ratio Tokamaks and is - selected via i_plasma_current=2. - - References: - - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). - 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), - 1729-1738. https://doi.org/10.13182/FST92-A29971 - """ - - # Transform q95 to qbar - qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 - - ff1, ff2, d1, d2 = _plascar_bpol(aspect, eps, kappa, delta) - - e1 = (2.0 * kappa) / (d1 * (1.0 + delta)) - e2 = (2.0 * kappa) / (d2 * (1.0 - delta)) - - return ( - rminor - * b_plasma_toroidal_on_axis - / qbar - * 5.0 - * kappa - / (2.0 * np.pi**2) - * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) - * (ff1 + ff2) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_ipdg89( - eps: float, kappa95: float, triang95: float -) -> float: - """ - Calculate the fq coefficient from the IPDG89 guidlines used in the plasma current scaling. - - Parameters: - - eps: float, plasma inverse aspect ratio - - kappa95: float, plasma elongation 95% - - triang95: float, plasma triangularity 95% - - Returns: - - float, the fq plasma current coefficient - - This function calculates the fq coefficient used in the IPDG89 plasma current scaling, - based on the given plasma parameters. - - References: - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - """ - return ( - 0.5 - * (1.17 - 0.65 * eps) - / ((1.0 - eps * eps) ** 2) - * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_todd( - eps: float, kappa95: float, triang95: float, model: int -) -> float: - """ - Calculate the fq coefficient used in the two Todd plasma current scalings. - - Parameters: - - eps: float, plasma inverse aspect ratio - - kappa95: float, plasma elongation 95% - - triang95: float, plasma triangularity 95% - - Returns: - - float, the fq plasma current coefficient - - This function calculates the fq coefficient based on the given plasma parameters for the two Todd scalings. - - References: - - D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - """ - # Calculate the Todd scaling based on the model - base_scaling = ( - (1.0 + 2.0 * eps**2) - * ((1.0 + kappa95**2) / 0.5) - * (1.24 - 0.54 * kappa95 + 0.3 * (kappa95**2 + triang95**2) + 0.125 * triang95) - ) - if model == 1: - return base_scaling - if model == 2: - return base_scaling * (1.0 + (abs(kappa95 - 1.2)) ** 3) - raise ProcessValueError(f"model = {model} is an invalid option") - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_hastie( - alphaj: float, - alphap: float, - b_plasma_toroidal_on_axis: float, - delta95: float, - eps: float, - kappa95: float, - pres_plasma_on_axis: float, - rmu0: float, -) -> float: - """ - Routine to calculate the f_q coefficient for the Connor-Hastie model used for scaling the plasma current. - - Parameters: - - alphaj: float, the current profile index - - alphap: float, the pressure profile index - - b_plasma_toroidal_on_axis: float, the toroidal field on axis (T) - - delta95: float, the plasma triangularity 95% - - eps: float, the inverse aspect ratio - - kappa95: float, the plasma elongation 95% - - pres_plasma_on_axis: float, the central plasma pressure (Pa) - - rmu0: float, the vacuum permeability (H/m) - - Returns: - - float, the F coefficient - - This routine calculates the f_q coefficient used for scaling the plasma current, - using the Connor-Hastie scaling - - Reference: - - J.W.Connor and R.J.Hastie, Culham Lab Report CLM-M106 (1985). - https://scientific-publications.ukaea.uk/wp-content/uploads/CLM-M106-1.pdf - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - """ - # Exponent in Connor-Hastie current profile - lamda = alphaj - - # Exponent in Connor-Hastie pressure profile - nu = alphap - - # Central plasma beta - beta0 = 2.0 * rmu0 * pres_plasma_on_axis / (b_plasma_toroidal_on_axis**2) - - # Plasma internal inductance - lamp1 = 1.0 + lamda - li = lamp1 / lamda * (lamp1 / lamda * np.log(lamp1) - 1.0) - - # T/r in AEA FUS 172 - kap1 = kappa95 + 1.0 - tr = kappa95 * delta95 / kap1**2 - - # E/r in AEA FUS 172 - er = (kappa95 - 1.0) / kap1 - - # T primed in AEA FUS 172 - tprime = 2.0 * tr * lamp1 / (1.0 + 0.5 * lamda) - - # E primed in AEA FUS 172 - eprime = er * lamp1 / (1.0 + lamda / 3.0) - - # Delta primed in AEA FUS 172 - deltap = (0.5 * kap1 * eps * 0.5 * li) + (beta0 / (0.5 * kap1 * eps)) * lamp1**2 / ( - 1.0 + nu - ) - - # Delta/R0 in AEA FUS 172 - deltar = beta0 / 6.0 * (1.0 + 5.0 * lamda / 6.0 + 0.25 * lamda**2) + ( - 0.5 * kap1 * eps - ) ** 2 * 0.125 * (1.0 - (lamda**2) / 3.0) - - # F coefficient - return (0.5 * kap1) ** 2 * ( - 1.0 - + eps**2 * (0.5 * kap1) ** 2 - + 0.5 * deltap**2 - + 2.0 * deltar - + 0.5 * (eprime**2 + er**2) - + 0.5 * (tprime**2 + 4.0 * tr**2) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_sauter( - eps: float, - kappa: float, - triang: float, -) -> float: - """ - Routine to calculate the f_q coefficient for the Sauter model used for scaling the plasma current. - - Parameters: - - eps: float, inverse aspect ratio - - kappa: float, plasma elongation at the separatrix - - triang: float, plasma triangularity at the separatrix - - Returns: - - float, the fq coefficient - - Reference: - - O. Sauter, Geometric formulas for system codes including the effect of negative triangularity, - Fusion Engineering and Design, Volume 112, 2016, Pages 633-645, - ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2016.04.033. - """ - w07 = 1.0 # zero squareness - can be modified later if required - - return ( - (4.1e6 / 5.0e6) - * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) - * (1.0 + 0.09 * triang + 0.16 * triang**2) - * (1.0 + 0.45 * triang * eps) - / (1.0 - 0.74 * eps) - * (1.0 + 0.55 * (w07 - 1.0)) - ) - - -@nb.jit(nopython=True, cache=True) -def calculate_current_coefficient_fiesta( - eps: float, kappa: float, triang: float -) -> float: - """ - Calculate the fq coefficient used in the FIESTA plasma current scaling. - - Parameters: - - eps: float, plasma inverse aspect ratio - - kappa: float, plasma elongation at the separatrix - - triang: float, plasma triangularity at the separatrix - - Returns: - - float, the fq plasma current coefficient - - This function calculates the fq coefficient based on the given plasma parameters for the FIESTA scaling. - - References: - - S.Muldrew et.al,“PROCESS”: Systems studies of spherical tokamaks, Fusion Engineering and Design, - Volume 154, 2020, 111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530. - """ - return 0.538 * (1.0 + 2.440 * eps**2.736) * kappa**2.154 * triang**0.060 - - # -------------------------------- # Bootstrap Current Calculations # -------------------------------- @@ -1457,7 +1170,7 @@ def physics(self): physics_variables.b_plasma_poloidal_average, physics_variables.qstar, physics_variables.plasma_current, - ) = self.calculate_plasma_current( + ) = self.current.calculate_plasma_current( physics_variables.alphaj, physics_variables.alphap, physics_variables.b_plasma_toroidal_on_axis, @@ -3291,238 +3004,72 @@ def plasma_ohmic_heating( return pden_plasma_ohmic_mw, p_plasma_ohmic_mw, f_res_plasma_neo, res_plasma - @staticmethod - def calculate_plasma_current( - alphaj: float, - alphap: float, - b_plasma_toroidal_on_axis: float, - eps: float, - i_plasma_current: int, - kappa: float, - kappa95: float, - pres_plasma_on_axis: float, - len_plasma_poloidal: float, - q95: float, - rmajor: float, - rminor: float, - triang: float, - triang95: float, - ) -> tuple[float, float, float, float, float]: - """Calculate the plasma current. - - Args: - alphaj (float): Current profile index. - alphap (float): Pressure profile index. - b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). - eps (float): Inverse aspect ratio. - i_plasma_current (int): Current scaling model to use. - 1 = Peng analytic fit - 2 = Peng divertor scaling (TART,STAR) - 3 = Simple ITER scaling - 4 = IPDG89 scaling - 5 = Todd empirical scaling I - 6 = Todd empirical scaling II - 7 = Connor-Hastie model - 8 = Sauter scaling (allowing negative triangularity) - 9 = FIESTA ST scaling - kappa (float): Plasma elongation. - kappa95 (float): Plasma elongation at 95% surface. - pres_plasma_on_axis (float): Central plasma pressure (Pa). - len_plasma_poloidal (float): Plasma perimeter length (m). - q95 (float): Plasma safety factor at 95% flux (= q-bar for i_plasma_current=2). - ind_plasma_internal_norm (float): Plasma normalised internal inductance. - rmajor (float): Major radius (m). - rminor (float): Minor radius (m). - triang (float): Plasma triangularity. - triang95 (float): Plasma triangularity at 95% surface. + def outtim(self): + po.oheadr(self.outfile, "Times") - Returns: - Tuple[float, float, float,]: Tuple containing b_plasma_poloidal_average, qstar, plasma_current, + po.ovarrf( + self.outfile, + "Initial charge time for CS from zero current (s)", + "(t_plant_pulse_coil_precharge)", + times_variables.t_plant_pulse_coil_precharge, + ) + po.ovarrf( + self.outfile, + "Plasma current ramp-up time (s)", + "(t_plant_pulse_plasma_current_ramp_up)", + times_variables.t_plant_pulse_plasma_current_ramp_up, + ) + po.ovarrf( + self.outfile, + "Heating time (s)", + "(t_plant_pulse_fusion_ramp)", + times_variables.t_plant_pulse_fusion_ramp, + ) + po.ovarre( + self.outfile, + "Burn time (s)", + "(t_plant_pulse_burn)", + times_variables.t_plant_pulse_burn, + "OP ", + ) + po.ovarrf( + self.outfile, + "Reset time to zero current for CS (s)", + "(t_plant_pulse_plasma_current_ramp_down)", + times_variables.t_plant_pulse_plasma_current_ramp_down, + ) + po.ovarrf( + self.outfile, + "Time between pulses (s)", + "(t_plant_pulse_dwell)", + times_variables.t_plant_pulse_dwell, + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Total plant cycle time (s)", + "(t_plant_pulse_total)", + times_variables.t_plant_pulse_total, + "OP ", + ) - Raises: - ValueError: If invalid value for i_plasma_current is provided. + def calculate_effective_charge_ionisation_profiles(self): + """Calculate the effective charge profiles for ionisation calculations.""" - Notes: - This routine calculates the plasma current based on the edge safety factor q95. - It will also make the current profile parameters consistent with the q-profile if required. - - References: - - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document - - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). - 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), - 1729-1738. https://doi.org/10.13182/FST92-A29971 - - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 - - M. Kovari et al, 2014, "PROCESS": A systems code for fusion power plants - Part 1: Physics - - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO - - T. Hartmann, 2013, Development of a modular systems code to analyse the implications of physics assumptions on the design of a demonstration fusion power plant - - Sauter, Geometric formulas for systems codes..., FED 2016 - """ - # Aspect ratio - aspect_ratio = 1.0 / eps - - # Only the Sauter scaling (i_plasma_current=8) is suitable for negative triangularity: - if i_plasma_current != 8 and triang < 0.0: - raise ProcessValueError( - f"Triangularity is negative without i_plasma_current = 8 selected: {triang=}, {i_plasma_current=}" - ) - - # Calculate the function Fq that scales the edge q from the - # circular cross-section cylindrical case - - # Peng analytical fit - if i_plasma_current == 1: - fq = calculate_current_coefficient_peng(eps, len_plasma_poloidal, rminor) - - # Peng scaling for double null divertor; TARTs [STAR Code] - elif i_plasma_current == 2: - plasma_current = 1.0e6 * calculate_plasma_current_peng( - q95, aspect_ratio, eps, rminor, b_plasma_toroidal_on_axis, kappa, triang - ) - - # Simple ITER scaling (simply the cylindrical case) - elif i_plasma_current == 3: - fq = 1.0 - - # ITER formula (IPDG89) - elif i_plasma_current == 4: - fq = calculate_current_coefficient_ipdg89(eps, kappa95, triang95) - - # Todd empirical scalings - # D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 - elif i_plasma_current in [5, 6]: - fq = calculate_current_coefficient_todd(eps, kappa95, triang95, model=1) - - if i_plasma_current == 6: - fq = calculate_current_coefficient_todd(eps, kappa95, triang95, model=2) - - # Connor-Hastie asymptotically-correct expression - elif i_plasma_current == 7: - fq = calculate_current_coefficient_hastie( - alphaj, - alphap, - b_plasma_toroidal_on_axis, - triang95, - eps, - kappa95, - pres_plasma_on_axis, - constants.RMU0, - ) - - # Sauter scaling allowing negative triangularity [FED May 2016] - # https://doi.org/10.1016/j.fusengdes.2016.04.033. - elif i_plasma_current == 8: - # Assumes zero squareness, note takes kappa, delta at separatrix not _95 - fq = calculate_current_coefficient_sauter(eps, kappa, triang) - - # FIESTA ST scaling - # https://doi.org/10.1016/j.fusengdes.2020.111530. - elif i_plasma_current == 9: - fq = calculate_current_coefficient_fiesta(eps, kappa, triang) - else: - raise ProcessValueError(f"Invalid value {i_plasma_current=}") - - # Main plasma current calculation using the fq value from the different settings - if i_plasma_current != 2: - plasma_current = ( - (2.0 * np.pi / constants.RMU0) - * rminor**2 - / (rmajor * q95) - * fq - * b_plasma_toroidal_on_axis - ) - # i_plasma_current == 2 case covered above - - # Calculate cyclindrical safety factor from IPDG89 - qstar = ( - 5.0e6 - * rminor**2 - / (rmajor * plasma_current / b_plasma_toroidal_on_axis) - * 0.5 - * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) - ) - - # Calculate the poloidal field generated by the plasma current - b_plasma_poloidal_average = calculate_poloidal_field( - i_plasma_current, - plasma_current, - q95, - aspect_ratio, - eps, - b_plasma_toroidal_on_axis, - kappa, - triang, - len_plasma_poloidal, - constants.RMU0, - ) - - return b_plasma_poloidal_average, qstar, plasma_current - - def outtim(self): - po.oheadr(self.outfile, "Times") - - po.ovarrf( - self.outfile, - "Initial charge time for CS from zero current (s)", - "(t_plant_pulse_coil_precharge)", - times_variables.t_plant_pulse_coil_precharge, - ) - po.ovarrf( - self.outfile, - "Plasma current ramp-up time (s)", - "(t_plant_pulse_plasma_current_ramp_up)", - times_variables.t_plant_pulse_plasma_current_ramp_up, - ) - po.ovarrf( - self.outfile, - "Heating time (s)", - "(t_plant_pulse_fusion_ramp)", - times_variables.t_plant_pulse_fusion_ramp, - ) - po.ovarre( - self.outfile, - "Burn time (s)", - "(t_plant_pulse_burn)", - times_variables.t_plant_pulse_burn, - "OP ", - ) - po.ovarrf( - self.outfile, - "Reset time to zero current for CS (s)", - "(t_plant_pulse_plasma_current_ramp_down)", - times_variables.t_plant_pulse_plasma_current_ramp_down, - ) - po.ovarrf( - self.outfile, - "Time between pulses (s)", - "(t_plant_pulse_dwell)", - times_variables.t_plant_pulse_dwell, - ) - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Total plant cycle time (s)", - "(t_plant_pulse_total)", - times_variables.t_plant_pulse_total, - "OP ", - ) - - def calculate_effective_charge_ionisation_profiles(self): - """Calculate the effective charge profiles for ionisation calculations.""" - - # Calculate the effective charge (zeff) profile across the plasma - # Returns an array of zeff at each radial point - zeff_profile = np.zeros_like(self.plasma_profile.teprofile.profile_y) - for i in range(len(zeff_profile)): - zeff_profile[i] = 0.0 - for imp in range(impurity_radiation_module.N_IMPURITIES): - zeff_profile[i] += ( - impurity_radiation_module.f_nd_impurity_electron_array[imp] - * impurity_radiation.zav_of_te( - imp, np.array([self.plasma_profile.teprofile.profile_y[i]]) - ).squeeze() - ** 2 - ) - physics_variables.n_charge_plasma_effective_profile = zeff_profile + # Calculate the effective charge (zeff) profile across the plasma + # Returns an array of zeff at each radial point + zeff_profile = np.zeros_like(self.plasma_profile.teprofile.profile_y) + for i in range(len(zeff_profile)): + zeff_profile[i] = 0.0 + for imp in range(impurity_radiation_module.N_IMPURITIES): + zeff_profile[i] += ( + impurity_radiation_module.f_nd_impurity_electron_array[imp] + * impurity_radiation.zav_of_te( + imp, np.array([self.plasma_profile.teprofile.profile_y[i]]) + ).squeeze() + ** 2 + ) + physics_variables.n_charge_plasma_effective_profile = zeff_profile # Assign the charge profiles of each species n_impurities = impurity_radiation_module.N_IMPURITIES @@ -9277,6 +8824,463 @@ def __init__(self): self.outfile = constants.NOUT self.mfile = constants.MFILE + def calculate_plasma_current( + self, + alphaj: float, + alphap: float, + b_plasma_toroidal_on_axis: float, + eps: float, + i_plasma_current: int, + kappa: float, + kappa95: float, + pres_plasma_on_axis: float, + len_plasma_poloidal: float, + q95: float, + rmajor: float, + rminor: float, + triang: float, + triang95: float, + ) -> tuple[float, float, float, float, float]: + """Calculate the plasma current. + + Args: + alphaj (float): Current profile index. + alphap (float): Pressure profile index. + b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). + eps (float): Inverse aspect ratio. + i_plasma_current (int): Current scaling model to use. + 1 = Peng analytic fit + 2 = Peng divertor scaling (TART,STAR) + 3 = Simple ITER scaling + 4 = IPDG89 scaling + 5 = Todd empirical scaling I + 6 = Todd empirical scaling II + 7 = Connor-Hastie model + 8 = Sauter scaling (allowing negative triangularity) + 9 = FIESTA ST scaling + kappa (float): Plasma elongation. + kappa95 (float): Plasma elongation at 95% surface. + pres_plasma_on_axis (float): Central plasma pressure (Pa). + len_plasma_poloidal (float): Plasma perimeter length (m). + q95 (float): Plasma safety factor at 95% flux (= q-bar for i_plasma_current=2). + ind_plasma_internal_norm (float): Plasma normalised internal inductance. + rmajor (float): Major radius (m). + rminor (float): Minor radius (m). + triang (float): Plasma triangularity. + triang95 (float): Plasma triangularity at 95% surface. + + Returns: + Tuple[float, float, float,]: Tuple containing b_plasma_poloidal_average, qstar, plasma_current, + + Raises: + ValueError: If invalid value for i_plasma_current is provided. + + Notes: + This routine calculates the plasma current based on the edge safety factor q95. + It will also make the current profile parameters consistent with the q-profile if required. + + References: + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729-1738. https://doi.org/10.13182/FST92-A29971 + - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + - M. Kovari et al, 2014, "PROCESS": A systems code for fusion power plants - Part 1: Physics + - H. Zohm et al, 2013, On the Physics Guidelines for a Tokamak DEMO + - T. Hartmann, 2013, Development of a modular systems code to analyse the implications of physics assumptions on the design of a demonstration fusion power plant + - Sauter, Geometric formulas for systems codes..., FED 2016 + """ + # Aspect ratio + aspect_ratio = 1.0 / eps + + # Only the Sauter scaling (i_plasma_current=8) is suitable for negative triangularity: + if i_plasma_current != 8 and triang < 0.0: + raise ProcessValueError( + f"Triangularity is negative without i_plasma_current = 8 selected: {triang=}, {i_plasma_current=}" + ) + + # Calculate the function Fq that scales the edge q from the + # circular cross-section cylindrical case + + # Peng analytical fit + if i_plasma_current == 1: + fq = self.current.calculate_current_coefficient_peng( + eps, len_plasma_poloidal, rminor + ) + + # Peng scaling for double null divertor; TARTs [STAR Code] + elif i_plasma_current == 2: + plasma_current = 1.0e6 * self.current.calculate_plasma_current_peng( + q95, aspect_ratio, eps, rminor, b_plasma_toroidal_on_axis, kappa, triang + ) + + # Simple ITER scaling (simply the cylindrical case) + elif i_plasma_current == 3: + fq = 1.0 + + # ITER formula (IPDG89) + elif i_plasma_current == 4: + fq = self.calculate_current_coefficient_ipdg89(eps, kappa95, triang95) + + # Todd empirical scalings + # D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 + elif i_plasma_current in [5, 6]: + fq = self.calculate_current_coefficient_todd(eps, kappa95, triang95, model=1) + + if i_plasma_current == 6: + fq = self.calculate_current_coefficient_todd( + eps, kappa95, triang95, model=2 + ) + + # Connor-Hastie asymptotically-correct expression + elif i_plasma_current == 7: + fq = self.calculate_current_coefficient_hastie( + alphaj, + alphap, + b_plasma_toroidal_on_axis, + triang95, + eps, + kappa95, + pres_plasma_on_axis, + constants.RMU0, + ) + + # Sauter scaling allowing negative triangularity [FED May 2016] + # https://doi.org/10.1016/j.fusengdes.2016.04.033. + elif i_plasma_current == 8: + # Assumes zero squareness, note takes kappa, delta at separatrix not _95 + fq = self.calculate_current_coefficient_sauter(eps, kappa, triang) + + # FIESTA ST scaling + # https://doi.org/10.1016/j.fusengdes.2020.111530. + elif i_plasma_current == 9: + fq = self.calculate_current_coefficient_fiesta(eps, kappa, triang) + else: + raise ProcessValueError(f"Invalid value {i_plasma_current=}") + + # Main plasma current calculation using the fq value from the different settings + if i_plasma_current != 2: + plasma_current = ( + (2.0 * np.pi / constants.RMU0) + * rminor**2 + / (rmajor * q95) + * fq + * b_plasma_toroidal_on_axis + ) + # i_plasma_current == 2 case covered above + + # Calculate cyclindrical safety factor from IPDG89 + qstar = ( + 5.0e6 + * rminor**2 + / (rmajor * plasma_current / b_plasma_toroidal_on_axis) + * 0.5 + * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) + ) + + # Calculate the poloidal field generated by the plasma current + b_plasma_poloidal_average = calculate_poloidal_field( + i_plasma_current, + plasma_current, + q95, + aspect_ratio, + eps, + b_plasma_toroidal_on_axis, + kappa, + triang, + len_plasma_poloidal, + constants.RMU0, + ) + + return b_plasma_poloidal_average, qstar, plasma_current + + @staticmethod + def calculate_current_coefficient_peng( + eps: float, len_plasma_poloidal: float, rminor: float + ) -> float: + """ + Calculate the plasma current scaling coefficient for the Peng scaling from the STAR code. + + :param eps: Plasma inverse aspect ratio. + :type eps: float + :param len_plasma_poloidal: Plasma poloidal perimeter length [m]. + :type len_plasma_poloidal: float + :param rminor: Plasma minor radius [m]. + :type rminor: float + + :return: The plasma current scaling coefficient. + :rtype: float + + :references: None + """ + + return ( + (1.22 - 0.68 * eps) + / ((1.0 - eps * eps) ** 2) + * (len_plasma_poloidal / (2.0 * np.pi * rminor)) ** 2 + ) + + @staticmethod + def calculate_plasma_current_peng( + q95: float, + aspect: float, + eps: float, + rminor: float, + b_plasma_toroidal_on_axis: float, + kappa: float, + delta: float, + ) -> float: + """ + Function to calculate plasma current (Peng scaling from the STAR code) + + Parameters: + - q95: float, 95% flux surface safety factor + - aspect: float, plasma aspect ratio + - eps: float, inverse aspect ratio + - rminor: float, plasma minor radius (m) + - b_plasma_toroidal_on_axis: float, toroidal field on axis (T) + - kappa: float, plasma elongation + - delta: float, plasma triangularity + + Returns: + - float, plasma current in MA + + This function calculates the plasma current in MA, + using a scaling from Peng, Galambos and Shipe (1992). + It is primarily used for Tight Aspect Ratio Tokamaks and is + selected via i_plasma_current=2. + + References: + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + - Peng, Y. K. M., Galambos, J. D., & Shipe, P. C. (1992). + 'Small Tokamaks for Fusion Technology Testing'. Fusion Technology, 21(3P2A), + 1729-1738. https://doi.org/10.13182/FST92-A29971 + """ + + # Transform q95 to qbar + qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 + + ff1, ff2, d1, d2 = _plascar_bpol(aspect, eps, kappa, delta) + + e1 = (2.0 * kappa) / (d1 * (1.0 + delta)) + e2 = (2.0 * kappa) / (d2 * (1.0 - delta)) + + return ( + rminor + * b_plasma_toroidal_on_axis + / qbar + * 5.0 + * kappa + / (2.0 * np.pi**2) + * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) + * (ff1 + ff2) + ) + + @staticmethod + def calculate_current_coefficient_ipdg89( + eps: float, kappa95: float, triang95: float + ) -> float: + """ + Calculate the fq coefficient from the IPDG89 guidlines used in the plasma current scaling. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa95: float, plasma elongation 95% + - triang95: float, plasma triangularity 95% + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient used in the IPDG89 plasma current scaling, + based on the given plasma parameters. + + References: + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + return ( + 0.5 + * (1.17 - 0.65 * eps) + / ((1.0 - eps * eps) ** 2) + * (1.0 + kappa95**2 * (1.0 + 2.0 * triang95**2 - 1.2 * triang95**3)) + ) + + @staticmethod + def calculate_current_coefficient_todd( + eps: float, kappa95: float, triang95: float, model: int + ) -> float: + """ + Calculate the fq coefficient used in the two Todd plasma current scalings. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa95: float, plasma elongation 95% + - triang95: float, plasma triangularity 95% + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient based on the given plasma parameters for the two Todd scalings. + + References: + - D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + # Calculate the Todd scaling based on the model + base_scaling = ( + (1.0 + 2.0 * eps**2) + * ((1.0 + kappa95**2) / 0.5) + * ( + 1.24 + - 0.54 * kappa95 + + 0.3 * (kappa95**2 + triang95**2) + + 0.125 * triang95 + ) + ) + if model == 1: + return base_scaling + if model == 2: + return base_scaling * (1.0 + (abs(kappa95 - 1.2)) ** 3) + raise ProcessValueError(f"model = {model} is an invalid option") + + @staticmethod + def calculate_current_coefficient_hastie( + alphaj: float, + alphap: float, + b_plasma_toroidal_on_axis: float, + delta95: float, + eps: float, + kappa95: float, + pres_plasma_on_axis: float, + rmu0: float, + ) -> float: + """ + Routine to calculate the f_q coefficient for the Connor-Hastie model used for scaling the plasma current. + + Parameters: + - alphaj: float, the current profile index + - alphap: float, the pressure profile index + - b_plasma_toroidal_on_axis: float, the toroidal field on axis (T) + - delta95: float, the plasma triangularity 95% + - eps: float, the inverse aspect ratio + - kappa95: float, the plasma elongation 95% + - pres_plasma_on_axis: float, the central plasma pressure (Pa) + - rmu0: float, the vacuum permeability (H/m) + + Returns: + - float, the F coefficient + + This routine calculates the f_q coefficient used for scaling the plasma current, + using the Connor-Hastie scaling + + Reference: + - J.W.Connor and R.J.Hastie, Culham Lab Report CLM-M106 (1985). + https://scientific-publications.ukaea.uk/wp-content/uploads/CLM-M106-1.pdf + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + # Exponent in Connor-Hastie current profile + lamda = alphaj + + # Exponent in Connor-Hastie pressure profile + nu = alphap + + # Central plasma beta + beta0 = 2.0 * rmu0 * pres_plasma_on_axis / (b_plasma_toroidal_on_axis**2) + + # Plasma internal inductance + lamp1 = 1.0 + lamda + li = lamp1 / lamda * (lamp1 / lamda * np.log(lamp1) - 1.0) + + # T/r in AEA FUS 172 + kap1 = kappa95 + 1.0 + tr = kappa95 * delta95 / kap1**2 + + # E/r in AEA FUS 172 + er = (kappa95 - 1.0) / kap1 + + # T primed in AEA FUS 172 + tprime = 2.0 * tr * lamp1 / (1.0 + 0.5 * lamda) + + # E primed in AEA FUS 172 + eprime = er * lamp1 / (1.0 + lamda / 3.0) + + # Delta primed in AEA FUS 172 + deltap = (0.5 * kap1 * eps * 0.5 * li) + ( + beta0 / (0.5 * kap1 * eps) + ) * lamp1**2 / (1.0 + nu) + + # Delta/R0 in AEA FUS 172 + deltar = beta0 / 6.0 * (1.0 + 5.0 * lamda / 6.0 + 0.25 * lamda**2) + ( + 0.5 * kap1 * eps + ) ** 2 * 0.125 * (1.0 - (lamda**2) / 3.0) + + # F coefficient + return (0.5 * kap1) ** 2 * ( + 1.0 + + eps**2 * (0.5 * kap1) ** 2 + + 0.5 * deltap**2 + + 2.0 * deltar + + 0.5 * (eprime**2 + er**2) + + 0.5 * (tprime**2 + 4.0 * tr**2) + ) + + @staticmethod + def calculate_current_coefficient_sauter( + eps: float, + kappa: float, + triang: float, + ) -> float: + """ + Routine to calculate the f_q coefficient for the Sauter model used for scaling the plasma current. + + Parameters: + - eps: float, inverse aspect ratio + - kappa: float, plasma elongation at the separatrix + - triang: float, plasma triangularity at the separatrix + + Returns: + - float, the fq coefficient + + Reference: + - O. Sauter, Geometric formulas for system codes including the effect of negative triangularity, + Fusion Engineering and Design, Volume 112, 2016, Pages 633-645, + ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2016.04.033. + """ + w07 = 1.0 # zero squareness - can be modified later if required + + return ( + (4.1e6 / 5.0e6) + * (1.0 + 1.2 * (kappa - 1.0) + 0.56 * (kappa - 1.0) ** 2) + * (1.0 + 0.09 * triang + 0.16 * triang**2) + * (1.0 + 0.45 * triang * eps) + / (1.0 - 0.74 * eps) + * (1.0 + 0.55 * (w07 - 1.0)) + ) + + @staticmethod + def calculate_current_coefficient_fiesta( + eps: float, kappa: float, triang: float + ) -> float: + """ + Calculate the fq coefficient used in the FIESTA plasma current scaling. + + Parameters: + - eps: float, plasma inverse aspect ratio + - kappa: float, plasma elongation at the separatrix + - triang: float, plasma triangularity at the separatrix + + Returns: + - float, the fq plasma current coefficient + + This function calculates the fq coefficient based on the given plasma parameters for the FIESTA scaling. + + References: + - S.Muldrew et.al,“PROCESS”: Systems studies of spherical tokamaks, Fusion Engineering and Design, + Volume 154, 2020, 111530, ISSN 0920-3796, https://doi.org/10.1016/j.fusengdes.2020.111530. + """ + return 0.538 * (1.0 + 2.440 * eps**2.736) * kappa**2.154 * triang**0.060 + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 9992a24a7..d257c7bf4 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -24,9 +24,8 @@ DetailedPhysics, Physics, PlasmaBeta, + PlasmaCurrent, PlasmaInductance, - calculate_current_coefficient_hastie, - calculate_plasma_current_peng, calculate_poloidal_field, diamagnetic_fraction_hender, diamagnetic_fraction_scene, @@ -56,6 +55,7 @@ def physics(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaCurrent(), ) @@ -1216,21 +1216,23 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): monkeypatch.setattr(physics_variables, "beta_total_vol_avg", plasmacurrentparam.beta) - b_plasma_poloidal_average, qstar, plasma_current = physics.calculate_plasma_current( - i_plasma_current=plasmacurrentparam.i_plasma_current, - alphaj=plasmacurrentparam.alphaj, - alphap=plasmacurrentparam.alphap, - b_plasma_toroidal_on_axis=plasmacurrentparam.b_plasma_toroidal_on_axis, - eps=plasmacurrentparam.eps, - kappa=plasmacurrentparam.kappa, - kappa95=plasmacurrentparam.kappa95, - pres_plasma_on_axis=plasmacurrentparam.pres_plasma_on_axis, - len_plasma_poloidal=plasmacurrentparam.len_plasma_poloidal, - q95=plasmacurrentparam.q95, - rmajor=plasmacurrentparam.rmajor, - rminor=plasmacurrentparam.rminor, - triang=plasmacurrentparam.triang, - triang95=plasmacurrentparam.triang95, + b_plasma_poloidal_average, qstar, plasma_current = ( + PlasmaCurrent().calculate_plasma_current( + i_plasma_current=plasmacurrentparam.i_plasma_current, + alphaj=plasmacurrentparam.alphaj, + alphap=plasmacurrentparam.alphap, + b_plasma_toroidal_on_axis=(plasmacurrentparam.b_plasma_toroidal_on_axis), + eps=plasmacurrentparam.eps, + kappa=plasmacurrentparam.kappa, + kappa95=plasmacurrentparam.kappa95, + pres_plasma_on_axis=plasmacurrentparam.pres_plasma_on_axis, + len_plasma_poloidal=plasmacurrentparam.len_plasma_poloidal, + q95=plasmacurrentparam.q95, + rmajor=plasmacurrentparam.rmajor, + rminor=plasmacurrentparam.rminor, + triang=plasmacurrentparam.triang, + triang95=plasmacurrentparam.triang95, + ) ) assert b_plasma_poloidal_average == pytest.approx(plasmacurrentparam.expected_bp) @@ -1270,7 +1272,9 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): ), ) def test_calculate_plasma_current_peng(arguments, expected): - assert calculate_plasma_current_peng(**arguments) == pytest.approx(expected) + assert PlasmaCurrent.calculate_plasma_current_peng(**arguments) == pytest.approx( + expected + ) @pytest.mark.parametrize( @@ -1334,7 +1338,7 @@ def test_calculate_beta_limit(): def test_conhas(): - assert calculate_current_coefficient_hastie( + assert PlasmaCurrent.calculate_current_coefficient_hastie( 5, 5, 12, 0.5, 0.33, 1.85, 2e3, constants.RMU0 ) == pytest.approx(2.518876726889116) diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index b72371607..6a2b44052 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -28,7 +28,7 @@ ) from process.fw import Fw from process.hcpb import CCFE_HCPB -from process.physics import Physics, PlasmaBeta, PlasmaInductance +from process.physics import Physics, PlasmaBeta, PlasmaCurrent, PlasmaInductance from process.plasma_profiles import PlasmaProfile from process.power import Power from process.stellarator import Neoclassics, Stellarator @@ -70,10 +70,12 @@ def stellarator(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaCurrent(), ), Neoclassics(), plasma_beta=PlasmaBeta(), plasma_inductance=PlasmaInductance(), + plasma_current=PlasmaCurrent(), )