From 1ddb045752dbfa490074af4234dda9698b4cdfc8 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 21:31:56 +0000 Subject: [PATCH 1/9] Add PlasmaCurrent class and integrate into Physics and Stellarator models --- process/main.py | 2 ++ process/models/physics/physics.py | 10 ++++++++++ 2 files changed, 12 insertions(+) diff --git a/process/main.py b/process/main.py index a2b243d10..a9a26ef12 100644 --- a/process/main.py +++ b/process/main.py @@ -706,6 +706,7 @@ def __init__(self): self.plasma_bootstrap_current = PlasmaBootstrapCurrent( plasma_profile=self.plasma_profile ) + self.plasma_current = PlasmaCurrent() self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive, @@ -714,6 +715,7 @@ def __init__(self): plasma_density_limit=self.plasma_density_limit, plasma_exhaust=self.plasma_exhaust, plasma_bootstrap_current=self.plasma_bootstrap_current, + plasma_current=self.plasma_current, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 84c686704..cf22c4c66 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -641,6 +641,7 @@ def __init__( plasma_density_limit, plasma_exhaust, plasma_bootstrap_current: PlasmaBootstrapCurrent, + plasma_current, ): self.outfile = constants.NOUT self.mfile = constants.MFILE @@ -651,6 +652,7 @@ def __init__( self.density_limit = plasma_density_limit self.exhaust = plasma_exhaust self.plasma_bootstrap_current = plasma_bootstrap_current + self.current = plasma_current def physics(self): """Routine to calculate tokamak plasma physics information @@ -7077,6 +7079,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.""" From afd75145e6c8f662310ad538898b63cbe864c666 Mon Sep 17 00:00:00 2001 From: chris-ashe Date: Sun, 8 Feb 2026 22:08:39 +0000 Subject: [PATCH 2/9] 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/models/physics/physics.py | 459 +++++++++++++++++++++++++++++- tests/unit/test_physics.py | 42 +-- tests/unit/test_stellarator.py | 1 + 3 files changed, 482 insertions(+), 20 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index cf22c4c66..745f3d908 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -706,7 +706,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, @@ -7086,6 +7086,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 52f5799ef..e479514fc 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -27,9 +27,8 @@ DetailedPhysics, Physics, PlasmaBeta, + PlasmaCurrent, PlasmaInductance, - calculate_current_coefficient_hastie, - calculate_plasma_current_peng, calculate_poloidal_field, diamagnetic_fraction_hender, diamagnetic_fraction_scene, @@ -62,6 +61,7 @@ def physics(): PlasmaDensityLimit(), PlasmaExhaust(), PlasmaBootstrapCurrent(plasma_profile=PlasmaProfile()), + PlasmaCurrent(), ) @@ -1224,21 +1224,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) @@ -1278,7 +1280,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( @@ -1342,7 +1346,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 6bec26a1f..3be36f256 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -92,6 +92,7 @@ def stellarator(): PlasmaDensityLimit(), PlasmaExhaust(), PlasmaBootstrapCurrent(plasma_profile=PlasmaProfile()), + PlasmaCurrent(), ), Neoclassics(), plasma_beta=PlasmaBeta(), From 76a67a5144b83ad5fe9697ebcab5339931185e94 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 26 Feb 2026 15:35:48 +0000 Subject: [PATCH 3/9] Create plasma current class --- process/main.py | 1 + process/models/physics/physics.py | 602 +--------------------- process/models/physics/plasma_current.py | 609 +++++++++++++++++++++++ tests/unit/test_physics.py | 7 +- tests/unit/test_stellarator.py | 1 + 5 files changed, 617 insertions(+), 603 deletions(-) create mode 100644 process/models/physics/plasma_current.py diff --git a/process/main.py b/process/main.py index a9a26ef12..9a8f6fc2f 100644 --- a/process/main.py +++ b/process/main.py @@ -102,6 +102,7 @@ PlasmaBeta, PlasmaInductance, ) +from process.models.physics.plasma_current import PlasmaCurrent from process.models.physics.plasma_geometry import PlasmaGeom from process.models.physics.plasma_profiles import PlasmaProfile from process.models.power import Power diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 745f3d908..c429db43c 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -89,140 +89,6 @@ def rether( # ----------------------------------------------------- -@nb.jit(nopython=True, cache=True) -def _plascar_bpol( - aspect: float, eps: float, kappa: float, delta: float -) -> tuple[float, float, float, float]: - """Calculate the poloidal field coefficients for determining the plasma current - and poloidal field. - - - This internal function calculates the poloidal field coefficients, - which is used to calculate the poloidal field and the plasma current. - - Parameters - ---------- - aspect : - plasma aspect ratio - eps : - inverse aspect ratio - kappa : - plasma elongation - delta : - plasma triangularity - - Returns - ------- - : - coefficients ff1, ff2, d1, d2 - - References - ---------- - - 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 - - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - - """ - # Original coding, only suitable for TARTs [STAR Code] - - c1 = (kappa**2 / (1.0 + delta)) + delta - c2 = (kappa**2 / (1.0 - delta)) - delta - - d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 - d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 - - c1_aspect = ((c1 * eps) - 1.0) if aspect < c1 else (1.0 - (c1 * eps)) - - y1 = np.sqrt(c1_aspect / (1.0 + eps)) * ((1.0 + delta) / kappa) - y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * ((1.0 - delta) / kappa) - - h2 = (1.0 + (c2 - 1.0) * (eps / 2.0)) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) - f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * ((c2 * eps) + 1.0)) - g = (eps * kappa) / (1.0 - (eps * delta)) - ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) - - h1 = (1.0 + (1.0 - c1) * (eps / 2.0)) / np.sqrt((1.0 + eps) * c1_aspect) - f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) - - if aspect < c1: - ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) - else: - ff1 = -f1 * (-g + 2.0 * h1 * np.arctan(y1)) - - return ff1, ff2, d1, d2 - - -def calculate_poloidal_field( - i_plasma_current: int, - ip: float, - q95: float, - aspect: float, - eps: float, - b_plasma_toroidal_on_axis: float, - kappa: float, - delta: float, - perim: float, - rmu0: float, -) -> float: - """Function to calculate poloidal field from the plasma current - - This function calculates the poloidal field from the plasma current in Tesla, - using a simple calculation using Ampere's law for conventional - tokamaks, or for TARTs, a scaling from Peng, Galambos and - Shipe (1992). - - Parameters - ---------- - i_plasma_current : - current scaling model to use - ip : - plasma current (A) - q95 : - 95% flux surface safety factor - aspect : - plasma aspect ratio - eps : - inverse aspect ratio - b_plasma_toroidal_on_axis : - toroidal field on axis (T) - kappa : - plasma elongation - delta : - plasma triangularity - perim : - plasma perimeter (m) - rmu0 : - vacuum permeability (H/m) - - Returns - ------- - : - poloidal field in Tesla - - - 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 - - """ - # Use Ampere's law using the plasma poloidal cross-section - if i_plasma_current != 2: - return rmu0 * ip / perim - # Use the relation from Peng, Galambos and Shipe (1992) [STAR code] otherwise - ff1, ff2, _, _ = _plascar_bpol(aspect, eps, kappa, delta) - - # Transform q95 to qbar - qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 - - 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: @@ -252,6 +118,7 @@ def calculate_current_coefficient_peng( def calculate_plasma_current_peng( + self, q95: float, aspect: float, eps: float, @@ -304,7 +171,7 @@ def calculate_plasma_current_peng( # 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) + ff1, ff2, d1, d2 = self._plascar_bpol(aspect, eps, kappa, delta) e1 = (2.0 * kappa) / (d1 * (1.0 + delta)) e2 = (2.0 * kappa) / (d2 * (1.0 - delta)) @@ -7079,471 +6946,6 @@ 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 - - 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/process/models/physics/plasma_current.py b/process/models/physics/plasma_current.py new file mode 100644 index 000000000..c8a3984aa --- /dev/null +++ b/process/models/physics/plasma_current.py @@ -0,0 +1,609 @@ +import logging + +import numpy as np + +from process import constants +from process.data_structure import ( + physics_variables, +) +from process.exceptions import ProcessValueError + +logger = logging.getLogger(__name__) + + +class PlasmaCurrent: + """Class to hold plasma current calculations for plasma processing.""" + + 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 = self.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 _plascar_bpol( + aspect: float, eps: float, kappa: float, delta: float + ) -> tuple[float, float, float, float]: + """Calculate the poloidal field coefficients for determining the plasma current + and poloidal field. + + + This internal function calculates the poloidal field coefficients, + which is used to calculate the poloidal field and the plasma current. + + Parameters + ---------- + aspect : + plasma aspect ratio + eps : + inverse aspect ratio + kappa : + plasma elongation + delta : + plasma triangularity + + Returns + ------- + : + coefficients ff1, ff2, d1, d2 + + References + ---------- + - 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 + - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + + """ + # Original coding, only suitable for TARTs [STAR Code] + + c1 = (kappa**2 / (1.0 + delta)) + delta + c2 = (kappa**2 / (1.0 - delta)) - delta + + d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 + d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 + + c1_aspect = ((c1 * eps) - 1.0) if aspect < c1 else (1.0 - (c1 * eps)) + + y1 = np.sqrt(c1_aspect / (1.0 + eps)) * ((1.0 + delta) / kappa) + y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * ((1.0 - delta) / kappa) + + h2 = (1.0 + (c2 - 1.0) * (eps / 2.0)) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) + f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * ((c2 * eps) + 1.0)) + g = (eps * kappa) / (1.0 - (eps * delta)) + ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) + + h1 = (1.0 + (1.0 - c1) * (eps / 2.0)) / np.sqrt((1.0 + eps) * c1_aspect) + f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) + + if aspect < c1: + ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) + else: + ff1 = -f1 * (-g + 2.0 * h1 * np.arctan(y1)) + + return ff1, ff2, d1, d2 + + def calculate_poloidal_field( + self, + i_plasma_current: int, + ip: float, + q95: float, + aspect: float, + eps: float, + b_plasma_toroidal_on_axis: float, + kappa: float, + delta: float, + perim: float, + rmu0: float, + ) -> float: + """Function to calculate poloidal field from the plasma current + + This function calculates the poloidal field from the plasma current in Tesla, + using a simple calculation using Ampere's law for conventional + tokamaks, or for TARTs, a scaling from Peng, Galambos and + Shipe (1992). + + Parameters + ---------- + i_plasma_current : + current scaling model to use + ip : + plasma current (A) + q95 : + 95% flux surface safety factor + aspect : + plasma aspect ratio + eps : + inverse aspect ratio + b_plasma_toroidal_on_axis : + toroidal field on axis (T) + kappa : + plasma elongation + delta : + plasma triangularity + perim : + plasma perimeter (m) + rmu0 : + vacuum permeability (H/m) + + Returns + ------- + : + poloidal field in Tesla + + + 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 + + """ + # Use Ampere's law using the plasma poloidal cross-section + if i_plasma_current != 2: + return rmu0 * ip / perim + # Use the relation from Peng, Galambos and Shipe (1992) [STAR code] otherwise + ff1, ff2, _, _ = self._plascar_bpol(aspect, eps, kappa, delta) + + # Transform q95 to qbar + qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 + + return b_plasma_toroidal_on_axis * (ff1 + ff2) / (2.0 * np.pi * qbar) + + @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 = PlasmaCurrent()._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 diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index e479514fc..6d4525abe 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -27,15 +27,14 @@ DetailedPhysics, Physics, PlasmaBeta, - PlasmaCurrent, PlasmaInductance, - calculate_poloidal_field, diamagnetic_fraction_hender, diamagnetic_fraction_scene, ps_fraction_scene, res_diff_time, rether, ) +from process.models.physics.plasma_current import PlasmaCurrent from process.models.physics.plasma_profiles import PlasmaProfile @@ -1336,7 +1335,9 @@ def test_calculate_plasma_current_peng(arguments, expected): ), ) def test_calculate_poloidal_field(arguments, expected): - assert calculate_poloidal_field(**arguments) == pytest.approx(expected) + assert PlasmaCurrent().calculate_poloidal_field(**arguments) == pytest.approx( + expected + ) def test_calculate_beta_limit(): diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 3be36f256..43ea5a0d1 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -37,6 +37,7 @@ PlasmaBeta, PlasmaInductance, ) +from process.models.physics.plasma_current import PlasmaCurrent from process.models.physics.plasma_profiles import PlasmaProfile from process.models.power import Power from process.models.stellarator.build import st_build From 632de115f014dd2658b6d16c4a48a36803050dd1 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 26 Feb 2026 15:51:10 +0000 Subject: [PATCH 4/9] Add cylindrical safety factor calculation to Physics model --- process/models/physics/physics.py | 55 +++++++++++++++++++++++++++++++ 1 file changed, 55 insertions(+) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index c429db43c..ac7c9b765 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -31,6 +31,52 @@ logger = logging.getLogger(__name__) +def calculate_cylindrical_safety_factor( + rmajor: float, + rminor: float, + plasma_current: float, + b_plasma_toroidal_on_axis: float, + kappa95: float, + triang95: float, +) -> float: + """Calculate the cylindrical safety factor from the IPDG89 guidelines. + + Parameters + ---------- + rmajor : float + Major radius of the tokamak in meters. + rminor : float + Minor radius of the tokamak in meters. + plasma_current : float + Plasma current in amperes. + b_plasma_toroidal_on_axis : float + Toroidal magnetic field on axis in tesla. + kappa95 : float + Elongation at 95% of the plasma boundary. + triang95 : float + Triangularity at 95% of the plasma boundary. + Returns + ------- + float + Cylindrical safety factor (dimensionless). + Notes + ----- + The cylindrical safety factor is calculated following the IPDG89 guidelines. + The formula accounts for plasma elongation and triangularity effects on the + safety factor through the kappa95 and triang95 parameters. + + """ + + # Calculate cyclindrical safety factor from IPDG89 + return ( + ((2 * np.pi) / constants.RMU0) + * 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)) + ) + + @nb.jit(nopython=True, cache=True) def rether( alphan, @@ -590,6 +636,15 @@ def physics(self): physics_variables.triang95, ) + physics_variables.qstar = calculate_cylindrical_safety_factor( + rmajor=physics_variables.rmajor, + rminor=physics_variables.rminor, + plasma_current=physics_variables.plasma_current, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + kappa95=physics_variables.kappa95, + triang95=physics_variables.triang95, + ) + # ----------------------------------------------------- # Plasma Current Profile # ----------------------------------------------------- From 333341a820b93ad16576321ac4c7f02a312f3efb Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 6 Mar 2026 14:46:04 +0000 Subject: [PATCH 5/9] Fix import statements in plasma_current.py for consistency --- process/models/physics/plasma_current.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/process/models/physics/plasma_current.py b/process/models/physics/plasma_current.py index c8a3984aa..80e713c6b 100644 --- a/process/models/physics/plasma_current.py +++ b/process/models/physics/plasma_current.py @@ -2,11 +2,11 @@ import numpy as np -from process import constants +from process.core import constants +from process.core.exceptions import ProcessValueError from process.data_structure import ( physics_variables, ) -from process.exceptions import ProcessValueError logger = logging.getLogger(__name__) From 2bce9218feb1e0f8914a54c8849aceb2ea53010e Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 6 Mar 2026 14:59:50 +0000 Subject: [PATCH 6/9] Add cylindrical plasma current calculation method to PlasmaCurrent class --- process/models/physics/plasma_current.py | 31 ++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/process/models/physics/plasma_current.py b/process/models/physics/plasma_current.py index 80e713c6b..a2793ca0e 100644 --- a/process/models/physics/plasma_current.py +++ b/process/models/physics/plasma_current.py @@ -188,6 +188,37 @@ def calculate_plasma_current( return b_plasma_poloidal_average, qstar, plasma_current + @staticmethod + def calculate_cyclindrical_plasma_current( + rminor: float, rmajor: float, q95: float, b_plasma_toroidal_on_axis: float + ) -> float: + """Calculate the plasma current for a cylindrical plasma. + + Parameters + ---------- + rminor : + plasma minor radius (m) + rmajor : + plasma major radius (m) + q95 : + plasma safety factor at 95% flux + b_plasma_toroidal_on_axis : + toroidal field on axis (T) + + Returns + ------- + float + plasma current (A) + + """ + + return ( + (2.0 * np.pi / constants.RMU0) + * rminor**2 + / (rmajor * q95) + * b_plasma_toroidal_on_axis + ) + @staticmethod def _plascar_bpol( aspect: float, eps: float, kappa: float, delta: float From b962626c3dd18b4132cc4f5142f8df84ee068b35 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 6 Mar 2026 15:06:49 +0000 Subject: [PATCH 7/9] Refactor plasma current calculation to use PlasmaCurrentModel enum for improved readability and maintainability --- process/models/physics/plasma_current.py | 129 ++++++++++++++--------- 1 file changed, 79 insertions(+), 50 deletions(-) diff --git a/process/models/physics/plasma_current.py b/process/models/physics/plasma_current.py index a2793ca0e..bf9dc7edd 100644 --- a/process/models/physics/plasma_current.py +++ b/process/models/physics/plasma_current.py @@ -1,4 +1,5 @@ import logging +from enum import IntEnum import numpy as np @@ -11,6 +12,18 @@ logger = logging.getLogger(__name__) +class PlasmaCurrentModel(IntEnum): + PENG_ANALYTIC_FIT = 1 + PENG_DIVERTOR_SCALING = 2 + ITER_SCALING = 3 + IPDG89_SCALING = 4 + TODD_EMPIRICAL_SCALING_I = 5 + TODD_EMPIRICAL_SCALING_II = 6 + CONNOR_HASTIE_MODEL = 7 + SAUTER_SCALING = 8 + FIESTA_ST_SCALING = 9 + + class PlasmaCurrent: """Class to hold plasma current calculations for plasma processing.""" @@ -92,65 +105,81 @@ def calculate_plasma_current( raise ProcessValueError( f"Triangularity is negative without i_plasma_current = 8 selected: {triang=}, {i_plasma_current=}" ) + try: + model = PlasmaCurrentModel(int(physics_variables.i_plasma_current)) + # Calculate the function Fq that scales the edge q from the + # circular cross-section cylindrical case + + # Peng analytical fit + if model == PlasmaCurrentModel.PENG_ANALYTIC_FIT: + fq = self.current.calculate_current_coefficient_peng( + eps, len_plasma_poloidal, rminor + ) - # 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 + # Peng scaling for double null divertor; TARTs [STAR Code] + elif model == PlasmaCurrentModel.PENG_DIVERTOR_SCALING: + plasma_current = 1.0e6 * self.current.calculate_plasma_current_peng( + q95, + aspect_ratio, + eps, + rminor, + b_plasma_toroidal_on_axis, + kappa, + triang, + ) - # ITER formula (IPDG89) - elif i_plasma_current == 4: - fq = self.calculate_current_coefficient_ipdg89(eps, kappa95, triang95) + # Simple ITER scaling (simply the cylindrical case) + elif model == PlasmaCurrentModel.ITER_SCALING: + fq = 1.0 - # 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) + # ITER formula (IPDG89) + elif model == PlasmaCurrentModel.IPDG89_SCALING: + fq = self.calculate_current_coefficient_ipdg89(eps, kappa95, triang95) - if i_plasma_current == 6: + # Todd empirical scalings + # D.C.Robinson and T.N.Todd, Plasma and Contr Fusion 28 (1986) 1181 + elif model in [ + PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_I, + PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_II, + ]: fq = self.calculate_current_coefficient_todd( - eps, kappa95, triang95, model=2 + eps, kappa95, triang95, model=1 ) - # 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, - ) + if model == PlasmaCurrentModel.TODD_EMPIRICAL_SCALING_II: + fq = self.calculate_current_coefficient_todd( + eps, kappa95, triang95, model=2 + ) + + # Connor-Hastie asymptotically-correct expression + elif model == PlasmaCurrentModel.CONNOR_HASTIE_MODEL: + 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) + # Sauter scaling allowing negative triangularity [FED May 2016] + # https://doi.org/10.1016/j.fusengdes.2016.04.033. + elif model == PlasmaCurrentModel.SAUTER_SCALING: + # 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=}") + # FIESTA ST scaling + # https://doi.org/10.1016/j.fusengdes.2020.111530. + elif model == PlasmaCurrentModel.FIESTA_ST_SCALING: + fq = self.calculate_current_coefficient_fiesta(eps, kappa, triang) + + except ValueError: + raise ProcessValueError( + "Illegal value of i_plasma_current", + i_plasma_current=physics_variables.i_plasma_current, + ) from None # Main plasma current calculation using the fq value from the different settings if i_plasma_current != 2: From 72048b3a0095ec82fa218986ddb68d4651d18cb3 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 6 Mar 2026 15:11:42 +0000 Subject: [PATCH 8/9] :sparkle: Add additional plasma current variables to physics_variables module --- process/data_structure/physics_variables.py | 43 +++++++++++++++++++++ 1 file changed, 43 insertions(+) diff --git a/process/data_structure/physics_variables.py b/process/data_structure/physics_variables.py index 72e4fd564..1cee819fe 100644 --- a/process/data_structure/physics_variables.py +++ b/process/data_structure/physics_variables.py @@ -961,6 +961,31 @@ plasma_current: float = None """plasma current (A)""" +c_plasma_peng_analytic: float = None +"""Peng analytic plasma current (A)""" + +c_plasma_peng_double_null: float = None +"""Peng double null divertor plasma current (A)""" + +c_plasma_cyclindrical: float = None +"""Cylindrical plasma current (A)""" + +c_plasma_ipdg89: float = None +"""ITER IPDG89 plasma current (A)""" + +c_plasma_todd_empirical_i: float = None +"""Todd empirical plasma current I (A)""" + +c_plasma_todd_empirical_ii: float = None +"""Todd empirical plasma current II (A)""" +c_plasma_connor_hastie: float = None +"""Connor-Hastie plasma current (A)""" + +c_plasma_sauter: float = None +"""Sauter plasma current (A)""" + +c_plasma_fiesta_st: float = None +"""FIESTA ST plasma current (A)""" p_plasma_neutron_mw: float = None """Neutron fusion power from just the plasma [MW]""" @@ -1611,6 +1636,15 @@ def init_physics_variables(): pflux_fw_rad_mw, \ pden_ion_electron_equilibration_mw, \ plasma_current, \ + c_plasma_peng_analytic, \ + c_plasma_peng_double_null, \ + c_plasma_cyclindrical, \ + c_plasma_ipdg89, \ + c_plasma_todd_empirical_i, \ + c_plasma_todd_empirical_ii, \ + c_plasma_connor_hastie, \ + c_plasma_sauter, \ + c_plasma_fiesta_st, \ p_plasma_neutron_mw, \ p_neutron_total_mw, \ pden_neutron_total_mw, \ @@ -1893,6 +1927,15 @@ def init_physics_variables(): pflux_fw_rad_mw = 0.0 pden_ion_electron_equilibration_mw = 0.0 plasma_current = 0.0 + c_plasma_peng_analytic = 0.0 + c_plasma_peng_double_null = 0.0 + c_plasma_cyclindrical = 0.0 + c_plasma_ipdg89 = 0.0 + c_plasma_todd_empirical_i = 0.0 + c_plasma_todd_empirical_ii = 0.0 + c_plasma_connor_hastie = 0.0 + c_plasma_sauter = 0.0 + c_plasma_fiesta_st = 0.0 p_plasma_neutron_mw = 0.0 p_neutron_total_mw = 0.0 pden_neutron_total_mw = 0.0 From 4b3c5abf6dd11a7c5e17dcebf9ef205c8001a38f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Fri, 6 Mar 2026 15:35:16 +0000 Subject: [PATCH 9/9] Refactor plasma current calculations and update unit tests for improved clarity and maintainability --- process/models/physics/physics.py | 578 +---------------------- process/models/physics/plasma_current.py | 40 +- tests/unit/test_physics.py | 56 +-- 3 files changed, 47 insertions(+), 627 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index ac7c9b765..d993d5afd 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -130,357 +130,6 @@ def rether( return conie * (temp_plasma_ion_vol_avg_kev - te) / (te**1.5) -# ----------------------------------------------------- -# Plasma Current & Poloidal Field Calculations -# ----------------------------------------------------- - - -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. - - Parameters - ---------- - eps : float - Plasma inverse aspect ratio. - len_plasma_poloidal : float - Plasma poloidal perimeter length [m]. - rminor : float - Plasma minor radius [m]. - - Returns - ------- - float - The plasma current scaling coefficient. - - """ - - 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( - self, - 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) - - 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. - - Parameters - ---------- - q95 : - 95% flux surface safety factor - aspect : - plasma aspect ratio - eps : - inverse aspect ratio - rminor : - plasma minor radius (m) - b_plasma_toroidal_on_axis : - toroidal field on axis (T) - kappa : - plasma elongation - delta : - plasma triangularity - - - Returns - ------- - : - plasma current in MA - - - 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 = self._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. - - This function calculates the fq coefficient used in the IPDG89 plasma current scaling, - based on the given plasma parameters. - - Parameters - ---------- - eps : - plasma inverse aspect ratio - kappa95 : - plasma elongation 95% - triang95 : - plasma triangularity 95% - - Returns - ------- - : - the fq plasma current coefficient - - - 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. - - This function calculates the fq coefficient based on the given plasma parameters for the two Todd scalings. - - Parameters - ---------- - eps : - plasma inverse aspect ratio - kappa95 : - plasma elongation 95% - triang95 : - plasma triangularity 95% - model: - - Returns - ------- - : - the fq plasma current coefficient - - 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. - - This routine calculates the f_q coefficient used for scaling the plasma current, - using the Connor-Hastie scaling - - - Parameters - ---------- - alphaj : - the current profile index - alphap : - the pressure profile index - b_plasma_toroidal_on_axis : - the toroidal field on axis (T) - delta95 : - the plasma triangularity 95% - eps : - the inverse aspect ratio - kappa95 : - the plasma elongation 95% - pres_plasma_on_axis : - the central plasma pressure (Pa) - rmu0 : - the vacuum permeability (H/m) - - Returns - ------- - : - the F coefficient - - - 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 : - inverse aspect ratio - kappa : - plasma elongation at the separatrix - triang : - plasma triangularity at the separatrix - - Returns - ------- - : - 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. - - This function calculates the fq coefficient based on the given plasma parameters for the FIESTA scaling. - - Parameters - ---------- - eps : - plasma inverse aspect ratio - kappa : - plasma elongation at the separatrix - triang : - plasma triangularity at the separatrix - - Returns - ------- - : - the fq plasma current coefficient - - - 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 - - # ----------------------------------------------------- # Diamagnetic and Pfirsch-Schlüter Current Calculations # ----------------------------------------------------- @@ -615,11 +264,7 @@ def physics(self): ) # Calculate plasma current - ( - physics_variables.b_plasma_poloidal_average, - physics_variables.qstar, - physics_variables.plasma_current, - ) = self.current.calculate_plasma_current( + physics_variables.plasma_current = self.current.calculate_plasma_current( physics_variables.alphaj, physics_variables.alphap, physics_variables.b_plasma_toroidal_on_axis, @@ -645,6 +290,22 @@ def physics(self): triang95=physics_variables.triang95, ) + # Calculate the poloidal field generated by the plasma current + physics_variables.b_plasma_poloidal_average = ( + self.current.calculate_poloidal_field( + i_plasma_current=physics_variables.i_plasma_current, + ip=physics_variables.plasma_current, + q95=physics_variables.q95, + aspect=physics_variables.aspect, + eps=physics_variables.eps, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + kappa=physics_variables.kappa, + delta=physics_variables.triang, + perim=physics_variables.len_plasma_poloidal, + rmu0=constants.RMU0, + ) + ) + # ----------------------------------------------------- # Plasma Current Profile # ----------------------------------------------------- @@ -1950,211 +1611,6 @@ 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. - - 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. - - Parameters - ---------- - 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. - ValueError - If invalid value for i_plasma_current is provided. - ValueError - If invalid value for i_plasma_current is provided. - - References - ---------- - - J D Galambos, STAR Code - Spherical Tokamak Analysis and Reactor Code, unpublished internal Oak Ridge document - - 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 - - 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 - - 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 - - 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") diff --git a/process/models/physics/plasma_current.py b/process/models/physics/plasma_current.py index bf9dc7edd..59d93a8b0 100644 --- a/process/models/physics/plasma_current.py +++ b/process/models/physics/plasma_current.py @@ -184,38 +184,16 @@ def calculate_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) + self.calculate_cyclindrical_plasma_current( + rminor=rminor, + rmajor=rmajor, + q95=q95, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + ) * 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 = self.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 + return plasma_current @staticmethod def calculate_cyclindrical_plasma_current( @@ -407,8 +385,8 @@ def calculate_current_coefficient_peng( * (len_plasma_poloidal / (2.0 * np.pi * rminor)) ** 2 ) - @staticmethod def calculate_plasma_current_peng( + self, q95: float, aspect: float, eps: float, @@ -448,7 +426,7 @@ def calculate_plasma_current_peng( # Transform q95 to qbar qbar = q95 * 1.3e0 * (1.0e0 - physics_variables.eps) ** 0.6e0 - ff1, ff2, d1, d2 = PlasmaCurrent()._plascar_bpol(aspect, eps, kappa, delta) + ff1, ff2, d1, d2 = self._plascar_bpol(aspect, eps, kappa, delta) e1 = (2.0 * kappa) / (d1 * (1.0 + delta)) e2 = (2.0 * kappa) / (d2 * (1.0 - delta)) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 6d4525abe..0c9f6b3ad 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -1152,10 +1152,6 @@ class PlasmaCurrentParam(NamedTuple): expected_normalised_total_beta: Any = None - expected_bp: Any = None - - expected_qstar: Any = None - expected_plasma_current: Any = None @@ -1180,8 +1176,6 @@ class PlasmaCurrentParam(NamedTuple): triang=0.5, triang95=0.33333333333333331, expected_normalised_total_beta=2.4784688886891844, - expected_bp=0.96008591022564971, - expected_qstar=2.900802902105021, expected_plasma_current=18398455.678867526, ), PlasmaCurrentParam( @@ -1202,8 +1196,6 @@ class PlasmaCurrentParam(NamedTuple): triang=0.5, triang95=0.33333333333333331, expected_normalised_total_beta=2.4784688886891844, - expected_bp=0.96008591022564971, - expected_qstar=2.900802902105021, expected_plasma_current=18398455.678867526, ), ), @@ -1223,29 +1215,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 = ( - 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, - ) + plasma_current = physics.current.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) - - assert qstar == pytest.approx(plasmacurrentparam.expected_qstar) - assert plasma_current == pytest.approx(plasmacurrentparam.expected_plasma_current) @@ -1278,8 +1264,8 @@ def test_calculate_plasma_current(plasmacurrentparam, monkeypatch, physics): ), ), ) -def test_calculate_plasma_current_peng(arguments, expected): - assert PlasmaCurrent.calculate_plasma_current_peng(**arguments) == pytest.approx( +def test_calculate_plasma_current_peng(arguments, expected, physics): + assert physics.current.calculate_plasma_current_peng(**arguments) == pytest.approx( expected ) @@ -1334,8 +1320,8 @@ def test_calculate_plasma_current_peng(arguments, expected): ), ), ) -def test_calculate_poloidal_field(arguments, expected): - assert PlasmaCurrent().calculate_poloidal_field(**arguments) == pytest.approx( +def test_calculate_poloidal_field(arguments, expected, physics): + assert physics.current.calculate_poloidal_field(**arguments) == pytest.approx( expected ) @@ -1346,8 +1332,8 @@ def test_calculate_beta_limit(): ) == pytest.approx(0.0297619) -def test_conhas(): - assert PlasmaCurrent.calculate_current_coefficient_hastie( +def test_conhas(physics): + assert physics.current.calculate_current_coefficient_hastie( 5, 5, 12, 0.5, 0.33, 1.85, 2e3, constants.RMU0 ) == pytest.approx(2.518876726889116)