From a9660855f49a90ff2b2ee701307ab6f5d7205bac Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 10:41:02 +0000 Subject: [PATCH 01/19] Add PlasmaDensityLimit class and integrate into physics and stellarator models --- process/core/io/plot_proc.py | 3 + process/main.py | 2 + process/models/physics/physics.py | 207 +++++++++++++++++++++++++++--- tests/unit/test_physics.py | 4 +- tests/unit/test_stellarator.py | 1 + 5 files changed, 197 insertions(+), 20 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index af8b0dd94..42788f4c8 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -252,6 +252,9 @@ def plot_plasma( axis.plot(pg.rs[0], pg.zs[0], color="black") axis.plot(pg.rs[1], pg.zs[1], color="black") + axis.plot(pg.rs[0] + 0.005, pg.zs[0] + 0.005, color="red", linestyle="dashed") + axis.plot(pg.rs[1] + 0.005, pg.zs[1] + 0.005, color="red", linestyle="dashed") + # Set triang_95 to stop plotting plasma past boundary # Assume IPDG scaling triang_95 = triang / 1.5 diff --git a/process/main.py b/process/main.py index 8165ede28..9b504a21b 100644 --- a/process/main.py +++ b/process/main.py @@ -698,11 +698,13 @@ def __init__(self): ) self.plasma_beta = PlasmaBeta() self.plasma_inductance = PlasmaInductance() + self.plasma_density_limit = PlasmaDensityLimit() self.physics = Physics( plasma_profile=self.plasma_profile, current_drive=self.current_drive, plasma_beta=self.plasma_beta, plasma_inductance=self.plasma_inductance, + plasma_density_limit=self.plasma_density_limit, ) self.physics_detailed = DetailedPhysics( plasma_profile=self.plasma_profile, diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index d2405c263..6c511b7ea 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1625,13 +1625,14 @@ def _trapped_particle_fraction_sauter( class Physics: - def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance): + def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance, plasma_density_limit): self.outfile = constants.NOUT self.mfile = constants.MFILE self.plasma_profile = plasma_profile self.current_drive = current_drive self.beta = plasma_beta self.inductance = plasma_inductance + self.density_limit = plasma_density_limit def physics(self): """Routine to calculate tokamak plasma physics information @@ -2412,24 +2413,7 @@ def physics(self): physics_variables.kappa95, ) - # Density limit - ( - physics_variables.nd_plasma_electron_max_array, - physics_variables.nd_plasma_electrons_max, - ) = self.calculate_density_limit( - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.i_density_limit, - physics_variables.p_plasma_separatrix_mw, - current_drive_variables.p_hcd_injected_total_mw, - physics_variables.plasma_current, - divertor_variables.prn1, - physics_variables.qstar, - physics_variables.q95, - physics_variables.rmajor, - physics_variables.rminor, - physics_variables.a_plasma_surface, - physics_variables.n_charge_plasma_effective_vol_avg, - ) + self.density_limit.run() # Calculate transport losses and energy confinement time using the # chosen scaling law @@ -9727,6 +9711,191 @@ def output_volt_second_information(self): po.oblnkl(self.outfile) +class PlasmaDensityLimit: + """Class to hold plasma density limit calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + def run(self): + # Density limit + ( + physics_variables.nd_plasma_electron_max_array, + physics_variables.nd_plasma_electrons_max, + ) = self.calculate_density_limit( + physics_variables.b_plasma_toroidal_on_axis, + physics_variables.i_density_limit, + physics_variables.p_plasma_separatrix_mw, + current_drive_variables.p_hcd_injected_total_mw, + physics_variables.plasma_current, + divertor_variables.prn1, + physics_variables.qstar, + physics_variables.q95, + physics_variables.rmajor, + physics_variables.rminor, + physics_variables.a_plasma_surface, + physics_variables.n_charge_plasma_effective_vol_avg, + ) + + @staticmethod + def calculate_density_limit( + b_plasma_toroidal_on_axis: float, + i_density_limit: int, + p_plasma_separatrix_mw: float, + p_hcd_injected_total_mw: float, + plasma_current: float, + prn1: float, + qcyl: float, + q95: float, + rmajor: float, + rminor: float, + a_plasma_surface: float, + zeff: float, + ) -> tuple[np.ndarray, float]: + """ + Calculate the density limit using various models. + + Args: + b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). + i_density_limit (int): Switch denoting which formula to enforce. + p_plasma_separatrix_mw (float): Power flowing to the edge plasma via charged particles (MW). + p_hcd_injected_total_mw (float): Power injected into the plasma (MW). + plasma_current (float): Plasma current (A). + prn1 (float): Edge density / average plasma density. + qcyl (float): Equivalent cylindrical safety factor (qstar). + q95 (float): Safety factor at 95% surface. + rmajor (float): Plasma major radius (m). + rminor (float): Plasma minor radius (m). + a_plasma_surface (float): Plasma surface area (m^2). + zeff (float): Plasma effective charge. + + Returns: + Tuple[np.ndarray, float]: A tuple containing: + - nd_plasma_electron_max_array (np.ndarray): Average plasma density limit using seven different models (m^-3). + - nd_plasma_electrons_max (float): Enforced average plasma density limit (m^-3). + + Raises: + ValueError: If i_density_limit is not between 1 and 7. + + Notes: + This routine calculates several different formulae for the density limit and enforces the one chosen by the user. + For i_density_limit = 1-5, 8, we scale the sepatrix density limit output by the ratio of the separatrix to volume averaged density + + References: + - AEA FUS 172: Physics Assessment for the European Reactor Study + + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989 + + - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” + vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. ‌ + """ + + if i_density_limit < 1 or i_density_limit > 7: + raise ProcessValueError( + "Illegal value for i_density_limit", i_density_limit=i_density_limit + ) + + nd_plasma_electron_max_array = np.empty((8,)) + + # Power per unit area crossing the plasma edge + # (excludes radiation and neutrons) + + p_perp = p_plasma_separatrix_mw / a_plasma_surface + + # Old ASDEX density limit formula + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[0] = ( + 1.54e20 + * p_perp**0.43 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.45 + ) / prn1 + + # Borrass density limit model for ITER (I) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # Borrass et al, ITER-TN-PH-9-6 (1989) + + nd_plasma_electron_max_array[1] = ( + 1.8e20 + * p_perp**0.53 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.22 + ) / prn1 + + # Borrass density limit model for ITER (II) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # This formula is (almost) identical to that in the original routine + # denlim (now deleted). + + nd_plasma_electron_max_array[2] = ( + 0.5e20 + * p_perp**0.57 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.09 + ) / prn1 + + # JET edge radiation density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # qcyl=qstar here, but literature is not clear. + + denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) + if denom <= 0.0: + if i_density_limit == 4: + logger.error( + f"qcyl < 4/3; nd_plasma_electron_max_array(4) set to zero; model 5 will be enforced instead. {denom=} {qcyl=}" + ) + i_density_limit = 5 + + nd_plasma_electron_max_array[3] = 0.0 + else: + nd_plasma_electron_max_array[3] = ( + 1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom) + ) / prn1 + + # JET simplified density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[4] = ( + 0.237e20 + * b_plasma_toroidal_on_axis + * np.sqrt(p_plasma_separatrix_mw) + / rmajor + ) / prn1 + + # Hugill-Murakami M.q limit + # qcyl=qstar here, which is okay according to the literature + + nd_plasma_electron_max_array[5] = ( + 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) + ) + + # Greenwald limit + + nd_plasma_electron_max_array[6] = ( + 1.0e14 * plasma_current / (np.pi * rminor * rminor) + ) + + nd_plasma_electron_max_array[7] = ( + 1.0e20 + * 0.506 + * (p_hcd_injected_total_mw**0.396 * (plasma_current / 1.0e6) ** 0.265) + / (q95**0.323) + ) / prn1 + + # Enforce the chosen density limit + + return nd_plasma_electron_max_array, nd_plasma_electron_max_array[ + i_density_limit - 1 + ] + + 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 956912aa7..3b5f7a5c9 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -25,6 +25,7 @@ Physics, PlasmaBeta, PlasmaInductance, + PlasmaDensityLimit, calculate_current_coefficient_hastie, calculate_plasma_current_peng, calculate_poloidal_field, @@ -56,6 +57,7 @@ def physics(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaDensityLimit(), ) @@ -2473,7 +2475,7 @@ def test_calculate_density_limit(calculatedensitylimitparam, physics): """ nd_plasma_electron_max_array, nd_plasma_electrons_max = ( - physics.calculate_density_limit( + PlasmaDensityLimit().calculate_density_limit( i_density_limit=calculatedensitylimitparam.i_density_limit, b_plasma_toroidal_on_axis=calculatedensitylimitparam.b_plasma_toroidal_on_axis, p_plasma_separatrix_mw=calculatedensitylimitparam.p_plasma_separatrix_mw, diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 268a94b67..1cd6fca15 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -82,6 +82,7 @@ def stellarator(): ), PlasmaBeta(), PlasmaInductance(), + PlasmaDensityLimit(), ), Neoclassics(), plasma_beta=PlasmaBeta(), From 018d407939fe401836db4aa845e0133e772ad584 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 10:51:27 +0000 Subject: [PATCH 02/19] Add method to calculate Greenwald density limit in PlasmaDensityLimit class --- process/models/physics/physics.py | 30 ++++++++++++++++++++++++++++-- 1 file changed, 28 insertions(+), 2 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 6c511b7ea..e12379b41 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9739,7 +9739,33 @@ def run(self): ) @staticmethod + def calculate_greenwald_density_limit(c_plasma: float, rminor: float) -> float: + """ + Calculate the Greenwald density limit (n_GW). + + :param c_plasma: Plasma current (A). + :type c_plasma: float + :param rminor: Plasma minor radius (m). + :type rminor: float + :return: The Greenwald density limit (m⁻³). + :rtype: float + + :notes: The Greenwald limit is typically applied to the line averaged electron density + + :references: + - M. Greenwald et al., “A new look at density limits in tokamaks,” + Nuclear Fusion, vol. 28, no. 12, pp. 2199-2207, Dec. 1988, + doi: https://doi.org/10.1088/0029-5515/28/12/009. + + - M. Greenwald, “Density limits in toroidal plasmas,” + Plasma Physics and Controlled Fusion, vol. 44, no. 8, pp. R27-R53, Jul. 2002, + doi: https://doi.org/10.1088/0741-3335/44/8/201. + """ + + return 1.0e14 * c_plasma / (np.pi * rminor**2) + def calculate_density_limit( + self, b_plasma_toroidal_on_axis: float, i_density_limit: int, p_plasma_separatrix_mw: float, @@ -9878,8 +9904,8 @@ def calculate_density_limit( # Greenwald limit - nd_plasma_electron_max_array[6] = ( - 1.0e14 * plasma_current / (np.pi * rminor * rminor) + nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( + c_plasma=plasma_current, rminor=rminor ) nd_plasma_electron_max_array[7] = ( From a3a80b18847c89694033c3eac3d3d14e8404faad Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 11:01:00 +0000 Subject: [PATCH 03/19] Add ASDEX new density limit calculation method to PlasmaDensityLimit class --- .../source/physics-models/plasma_density.md | 4 +- process/models/physics/physics.py | 48 ++++++++++++++++--- 2 files changed, 44 insertions(+), 8 deletions(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index 5c4fa4ab3..ce44525cc 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -84,7 +84,7 @@ $$ ----------------- -## Greenwald model +## Greenwald model | `calculate_greenwald_density_limit()` Switch value: `i_density_limit = 7` [^3][^4] @@ -96,7 +96,7 @@ For the Greenwald model the limit applies to the line-averaged electron density, --------------------- -## ASDEX New model +## ASDEX New model | `calculate_asdex_new_density_limit()` Switch value: `i_density_limit = 8` [^5][^6] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index e12379b41..677cd9a19 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9764,6 +9764,42 @@ def calculate_greenwald_density_limit(c_plasma: float, rminor: float) -> float: return 1.0e14 * c_plasma / (np.pi * rminor**2) + @staticmethod + def calculate_asdex_new_density_limit( + p_hcd_injected_total_mw: float, c_plasma: float, q95: float, prn1: float + ) -> float: + """ + Calculate the ASDEX Upgrade new density limit. + + :param p_hcd_injected_total_mw: Power injected into the plasma (MW). + :type p_hcd_injected_total_mw: float + :param plasma_current: Plasma current (A). + :type plasma_current: float + :param q95: Safety factor at 95% surface. + :type q95: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The ASDEX Upgrade new density limit (m⁻³). + :rtype: float + + :notes: This limit is for the separatrix density so wee scale by `prn1` to get it as a volume average + + :references: + + - J. W. Berkery et al., “Density limits as disruption forecasters for spherical tokamaks,” + Plasma Physics and Controlled Fusion, vol. 65, no. 9, pp. 095003-095003, Jul. 2023, + doi: https://doi.org/10.1088/1361-6587/ace476. + + - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” vol. 57, no. 1, pp. 014038-014038, Nov. 2014, + doi: https://doi.org/10.1088/0741-3335/57/1/014038. + """ + return ( + 1.0e20 + * 0.506 + * (p_hcd_injected_total_mw**0.396 * (c_plasma / 1.0e6) ** 0.265) + / (q95**0.323) + ) / prn1 + def calculate_density_limit( self, b_plasma_toroidal_on_axis: float, @@ -9908,12 +9944,12 @@ def calculate_density_limit( c_plasma=plasma_current, rminor=rminor ) - nd_plasma_electron_max_array[7] = ( - 1.0e20 - * 0.506 - * (p_hcd_injected_total_mw**0.396 * (plasma_current / 1.0e6) ** 0.265) - / (q95**0.323) - ) / prn1 + nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( + p_hcd_injected_total_mw=p_hcd_injected_total_mw, + c_plasma=plasma_current, + q95=q95, + prn1=prn1, + ) # Enforce the chosen density limit From 5ebe78892d6e81010a01e9662f15afcc2b7ddbb2 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 11:08:52 +0000 Subject: [PATCH 04/19] Add Hugill-Murakami density limit calculation method to PlasmaDensityLimit class --- .../source/physics-models/plasma_density.md | 2 +- process/models/physics/physics.py | 26 +++++++++++++++++-- 2 files changed, 25 insertions(+), 3 deletions(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index ce44525cc..9c60d4d84 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -73,7 +73,7 @@ where $\kappa \approx 1.5, \Delta \approx 0.1a$ has been taken from JET. ----------------- -## Hugill-Murakami model +## Hugill-Murakami model | `calculate_hugill_murakami_density_limit()` Switch value: `i_density_limit = 6` [^2] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 677cd9a19..77c344226 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9738,6 +9738,28 @@ def run(self): physics_variables.n_charge_plasma_effective_vol_avg, ) + @staticmethod + def calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis: float, rmajor: float, qcyl: float + ) -> float: + """ + Calculate the Hugill-Murakami density limit. + + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param qcyl: Equivalent cylindrical safety factor (qstar). + :type qcyl: float + :return: The Hugill-Murakami density limit (m⁻³). + :rtype: float + + :references: + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', + """ + + return 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) + @staticmethod def calculate_greenwald_density_limit(c_plasma: float, rminor: float) -> float: """ @@ -9934,8 +9956,8 @@ def calculate_density_limit( # Hugill-Murakami M.q limit # qcyl=qstar here, which is okay according to the literature - nd_plasma_electron_max_array[5] = ( - 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) + nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, rmajor=rmajor, qcyl=qcyl ) # Greenwald limit From 2ceef94ebe21b689d904aa57fab7ce7f81c53ae5 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 11:19:08 +0000 Subject: [PATCH 05/19] Add JET simple density limit calculation method to PlasmaDensityLimit class --- .../source/physics-models/plasma_density.md | 2 +- process/models/physics/physics.py | 44 ++++++++++++++++--- 2 files changed, 39 insertions(+), 7 deletions(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index 9c60d4d84..cfceed67f 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -55,7 +55,7 @@ $$ ----------------- -## JET simplified model +## JET simplified model | `calculate_jet_simple_density_limit()` Switch value: `i_density_limit = 5` [^1] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 77c344226..503a07bf0 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9738,6 +9738,38 @@ def run(self): physics_variables.n_charge_plasma_effective_vol_avg, ) + @staticmethod + def calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis: float, + p_plasma_separatrix_mw: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the JET simple density limit. + + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param p_plasma_separatrix_mw: Power crossing the separatrix (MW). + :type p_plasma_separatrix_mw: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The JET simple density limit (m⁻³). + :rtype: float + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 0.237e20 + * b_plasma_toroidal_on_axis + * np.sqrt(p_plasma_separatrix_mw) + / rmajor + ) / prn1 + @staticmethod def calculate_hugill_murakami_density_limit( b_plasma_toroidal_on_axis: float, rmajor: float, qcyl: float @@ -9946,12 +9978,12 @@ def calculate_density_limit( # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. - nd_plasma_electron_max_array[4] = ( - 0.237e20 - * b_plasma_toroidal_on_axis - * np.sqrt(p_plasma_separatrix_mw) - / rmajor - ) / prn1 + nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + p_plasma_separatrix_mw=p_plasma_separatrix_mw, + rmajor=rmajor, + prn1=prn1, + ) # Hugill-Murakami M.q limit # qcyl=qstar here, which is okay according to the literature From eded56c1b21d5acb1011e68b88f00e8398a145e5 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 11:33:03 +0000 Subject: [PATCH 06/19] Add JET edge radiation density limit calculation method to PlasmaDensityLimit class --- .../source/physics-models/plasma_density.md | 2 +- process/models/physics/physics.py | 48 ++++++++++++++----- 2 files changed, 36 insertions(+), 14 deletions(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index cfceed67f..c6c2c1c84 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -45,7 +45,7 @@ $$ ----------------- -## JET edge radiation model +## JET edge radiation model | `calculate_jet_edge_radiation_density_limit()` Switch value: `i_density_limit = 4` [^1] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 503a07bf0..9cb40cab3 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9738,6 +9738,33 @@ def run(self): physics_variables.n_charge_plasma_effective_vol_avg, ) + @staticmethod + def calculate_jet_edge_radiation_density_limit( + zeff: float, p_hcd_injected_total_mw: float, prn1: float, qcyl: float + ) -> float: + """ + Calculate the JET edge radiation density limit. + + :param zeff: Effective charge (Z_eff). + :type zeff: float + :param p_hcd_injected_total_mw: Power injected into the plasma (MW). + :type p_hcd_injected_total_mw: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :param qcyl: Equivalent cylindrical safety factor (qstar). + :type qcyl: float + :return: The JET edge radiation density limit (m⁻³). + :rtype: float + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + + denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) + if denom <= 0.0: + return 0.0 + return (1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom)) / prn1 + @staticmethod def calculate_jet_simple_density_limit( b_plasma_toroidal_on_axis: float, @@ -9960,19 +9987,14 @@ def calculate_density_limit( # to give the density limit applying to the average plasma density. # qcyl=qstar here, but literature is not clear. - denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) - if denom <= 0.0: - if i_density_limit == 4: - logger.error( - f"qcyl < 4/3; nd_plasma_electron_max_array(4) set to zero; model 5 will be enforced instead. {denom=} {qcyl=}" - ) - i_density_limit = 5 - - nd_plasma_electron_max_array[3] = 0.0 - else: - nd_plasma_electron_max_array[3] = ( - 1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom) - ) / prn1 + nd_plasma_electron_max_array[3] = ( + self.calculate_jet_edge_radiation_density_limit( + zeff=zeff, + p_hcd_injected_total_mw=p_hcd_injected_total_mw, + prn1=prn1, + qcyl=qcyl, + ) + ) # JET simplified density limit model # This applies to the density at the plasma edge, so must be scaled From 97f925a0e4c703b7d0056f8566ce6b99e8d8991f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 11:38:13 +0000 Subject: [PATCH 07/19] Add Borrass ITER II density limit calculation method to PlasmaDensityLimit class --- .../source/physics-models/plasma_density.md | 2 +- process/models/physics/physics.py | 49 ++++++++++++++++--- 2 files changed, 44 insertions(+), 7 deletions(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index c6c2c1c84..aa9229ca1 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -35,7 +35,7 @@ $C \approx$ 1.8 for ITER-like conditions. ----------------- -## Borrass model for ITER, II +## Borrass model for ITER, II | `calculate_borrass_iter_ii_density_limit()` Switch value: `i_density_limit = 3` [^1] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 9cb40cab3..bd7e47fa8 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9738,6 +9738,42 @@ def run(self): physics_variables.n_charge_plasma_effective_vol_avg, ) + @staticmethod + def calculate_borrass_iter_ii_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the Borrass ITER II density limit. + + :param p_perp: Perpendicular power density (MW/m²). + :type p_perp: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param q95: Safety factor at 95% of the plasma poloidal flux. + :type q95: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The Borrass ITER II density limit (m⁻³). + :rtype: float + + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 0.5e20 + * p_perp**0.57 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.09 + ) / prn1 + @staticmethod def calculate_jet_edge_radiation_density_limit( zeff: float, p_hcd_injected_total_mw: float, prn1: float, qcyl: float @@ -9975,12 +10011,13 @@ def calculate_density_limit( # This formula is (almost) identical to that in the original routine # denlim (now deleted). - nd_plasma_electron_max_array[2] = ( - 0.5e20 - * p_perp**0.57 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.09 - ) / prn1 + nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) # JET edge radiation density limit model # This applies to the density at the plasma edge, so must be scaled From 05dc3e9693c76396ccc0138e3efa78efb26558f9 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 11:40:33 +0000 Subject: [PATCH 08/19] Add Borrass ITER I density limit calculation method to PlasmaDensityLimit class --- .../source/physics-models/plasma_density.md | 2 +- process/models/physics/physics.py | 49 ++++++++++++++++--- 2 files changed, 44 insertions(+), 7 deletions(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index aa9229ca1..0d3fa6272 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -23,7 +23,7 @@ $$ ----------------- -## Borrass model for ITER, I +## Borrass model for ITER, I | `calculate_borrass_iter_i_density_limit()` Switch value: `i_density_limit = 2` [^1] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index bd7e47fa8..75437fb02 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9738,6 +9738,42 @@ def run(self): physics_variables.n_charge_plasma_effective_vol_avg, ) + @staticmethod + def calculate_borrass_iter_i_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the Borrass ITER I density limit. + + :param p_perp: Perpendicular power density (MW/m²). + :type p_perp: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param q95: Safety factor at 95% of the plasma poloidal flux. + :type q95: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The Borrass ITER I density limit (m⁻³). + :rtype: float + + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 1.8e20 + * p_perp**0.53 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.22 + ) / prn1 + @staticmethod def calculate_borrass_iter_ii_density_limit( p_perp: float, @@ -9998,12 +10034,13 @@ def calculate_density_limit( # to give the density limit applying to the average plasma density. # Borrass et al, ITER-TN-PH-9-6 (1989) - nd_plasma_electron_max_array[1] = ( - 1.8e20 - * p_perp**0.53 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.22 - ) / prn1 + nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) # Borrass density limit model for ITER (II) # This applies to the density at the plasma edge, so must be scaled From d78da35f33b43f65232d394fd0db40dfa9c5de5f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 13:09:59 +0000 Subject: [PATCH 09/19] Add ASDEX density limit calculation method to PlasmaDensityLimit class --- .../source/physics-models/plasma_density.md | 2 +- process/models/physics/physics.py | 48 ++++++++++++++++--- 2 files changed, 43 insertions(+), 7 deletions(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index 0d3fa6272..04f6ca812 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -13,7 +13,7 @@ For the models below $P_{\perp}$ is the mean heat flux density across the separa ----------------- -## ASDEX model +## ASDEX model | `calculate_asdex_density_limit()` Switch value: `i_density_limit = 1`[^1][^2] diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 75437fb02..4c8e6be6e 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9738,6 +9738,41 @@ def run(self): physics_variables.n_charge_plasma_effective_vol_avg, ) + @staticmethod + def calculate_asdex_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the ASDEX density limit. + + :param p_perp: Perpendicular power density (MW/m²). + :type p_perp: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param q95: Safety factor at 95% of the plasma poloidal flux. + :type q95: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The ASDEX density limit (m⁻³). + :rtype: float + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 1.54e20 + * p_perp**0.43 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.45 + ) / prn1 + @staticmethod def calculate_borrass_iter_i_density_limit( p_perp: float, @@ -10022,12 +10057,13 @@ def calculate_density_limit( # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. - nd_plasma_electron_max_array[0] = ( - 1.54e20 - * p_perp**0.43 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.45 - ) / prn1 + nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) # Borrass density limit model for ITER (I) # This applies to the density at the plasma edge, so must be scaled From 3d66f17d5fea86a10320c7a1c07a1d230de5825c Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 13:23:18 +0000 Subject: [PATCH 10/19] Add DensityLimitModel class for electron density model types --- process/models/physics/physics.py | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 4c8e6be6e..b45ed7dfb 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9711,6 +9711,19 @@ def output_volt_second_information(self): po.oblnkl(self.outfile) +class DensityLimitModel(IntEnum): + """Electron density model types""" + + ASDEX = 1 + BORRASS_ITER_I = 2 + BORRASS_ITER_II = 3 + JET_EDGE_RADIATION = 4 + JET_SIMPLE = 5 + HUGILL_MURAKAMI = 6 + GREENWALD = 7 + ASDEX_NEW = 8 + + class PlasmaDensityLimit: """Class to hold plasma density limit calculations for plasma processing.""" From 60a1cbcd5ccd19c5c07b692c8c2c3620c4cb1ec0 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 13:41:00 +0000 Subject: [PATCH 11/19] Implement multiple density limit models in PlasmaDensityLimit class --- process/models/physics/physics.py | 148 ++++++++++++++++++++++++++---- 1 file changed, 131 insertions(+), 17 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index b45ed7dfb..70e5405c5 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -9732,25 +9732,139 @@ def __init__(self): self.mfile = constants.MFILE def run(self): - # Density limit - ( - physics_variables.nd_plasma_electron_max_array, - physics_variables.nd_plasma_electrons_max, - ) = self.calculate_density_limit( - physics_variables.b_plasma_toroidal_on_axis, - physics_variables.i_density_limit, - physics_variables.p_plasma_separatrix_mw, - current_drive_variables.p_hcd_injected_total_mw, - physics_variables.plasma_current, - divertor_variables.prn1, - physics_variables.qstar, - physics_variables.q95, - physics_variables.rmajor, - physics_variables.rminor, - physics_variables.a_plasma_surface, - physics_variables.n_charge_plasma_effective_vol_avg, + nd_plasma_electron_max_array = np.empty((8,)) + + p_perp = ( + physics_variables.p_plasma_separatrix_mw / physics_variables.a_plasma_surface ) + # Old ASDEX density limit formula + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # Borrass density limit model for ITER (I) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # Borrass et al, ITER-TN-PH-9-6 (1989) + + nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # Borrass density limit model for ITER (II) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # This formula is (almost) identical to that in the original routine + # denlim (now deleted). + + nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # JET edge radiation density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # qcyl=qstar here, but literature is not clear. + + nd_plasma_electron_max_array[3] = ( + self.calculate_jet_edge_radiation_density_limit( + zeff=physics_variables.n_charge_plasma_effective_vol_avg, + p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, + prn1=divertor_variables.prn1, + qcyl=physics_variables.qstar, + ) + ) + + # JET simplified density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + p_plasma_separatrix_mw=physics_variables.p_plasma_separatrix_mw, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # Hugill-Murakami M.q limit + # qcyl=qstar here, which is okay according to the literature + + nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + rmajor=physics_variables.rmajor, + qcyl=physics_variables.qstar, + ) + + # Greenwald limit + + nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( + c_plasma=physics_variables.plasma_current, rminor=physics_variables.rminor + ) + + nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( + p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, + c_plasma=physics_variables.plasma_current, + q95=physics_variables.q95, + prn1=divertor_variables.prn1, + ) + + physics_variables.nd_plasma_electron_max_array = nd_plasma_electron_max_array + + # Calculate beta_norm_max based on i_beta_norm_max + try: + model = DensityLimitModel(int(physics_variables.i_density_limit)) + physics_variables.nd_plasma_electrons_max = self.get_density_limit_value( + model + ) + except ValueError: + raise ProcessValueError( + "Illegal value of i_density_limit", + i_density_limit=physics_variables.i_density_limit, + ) from None + + def get_density_limit_value(self, model: DensityLimitModel) -> float: + """Get the density limit value (n_e_max) for the specified model.""" + model_map = { + DensityLimitModel.ASDEX: physics_variables.nd_plasma_electron_max_array[0], + DensityLimitModel.BORRASS_ITER_I: physics_variables.nd_plasma_electron_max_array[ + 1 + ], + DensityLimitModel.BORRASS_ITER_II: physics_variables.nd_plasma_electron_max_array[ + 2 + ], + DensityLimitModel.JET_EDGE_RADIATION: physics_variables.nd_plasma_electron_max_array[ + 3 + ], + DensityLimitModel.JET_SIMPLE: physics_variables.nd_plasma_electron_max_array[ + 4 + ], + DensityLimitModel.HUGILL_MURAKAMI: physics_variables.nd_plasma_electron_max_array[ + 5 + ], + DensityLimitModel.GREENWALD: physics_variables.nd_plasma_electron_max_array[ + 6 + ], + DensityLimitModel.ASDEX_NEW: physics_variables.nd_plasma_electron_max_array[ + 7 + ], + } + return model_map[model] + @staticmethod def calculate_asdex_density_limit( p_perp: float, From 819b7a9464529aad0b12d23a9c5e157f00c81384 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 12 Feb 2026 14:10:41 +0000 Subject: [PATCH 12/19] Refactor density limit output handling in Physics class --- process/models/physics/physics.py | 150 ++++++++++++++++-------------- 1 file changed, 81 insertions(+), 69 deletions(-) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index 70e5405c5..db7ff7c20 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -1625,7 +1625,14 @@ def _trapped_particle_fraction_sauter( class Physics: - def __init__(self, plasma_profile, current_drive, plasma_beta, plasma_inductance, plasma_density_limit): + def __init__( + self, + plasma_profile, + current_drive, + plasma_beta, + plasma_inductance, + plasma_density_limit, + ): self.outfile = constants.NOUT self.mfile = constants.MFILE self.plasma_profile = plasma_profile @@ -4570,74 +4577,7 @@ def outplas(self): po.oblnkl(self.outfile) if stellarator_variables.istell == 0: - po.osubhd(self.outfile, "Density Limit using different models :") - po.ovarre( - self.outfile, - "Old ASDEX model", - "(nd_plasma_electron_max_array(1))", - physics_variables.nd_plasma_electron_max_array[0], - "OP ", - ) - po.ovarre( - self.outfile, - "Borrass ITER model I", - "(nd_plasma_electron_max_array(2))", - physics_variables.nd_plasma_electron_max_array[1], - "OP ", - ) - po.ovarre( - self.outfile, - "Borrass ITER model II", - "(nd_plasma_electron_max_array(3))", - physics_variables.nd_plasma_electron_max_array[2], - "OP ", - ) - po.ovarre( - self.outfile, - "JET edge radiation model", - "(nd_plasma_electron_max_array(4))", - physics_variables.nd_plasma_electron_max_array[3], - "OP ", - ) - po.ovarre( - self.outfile, - "JET simplified model", - "(nd_plasma_electron_max_array(5))", - physics_variables.nd_plasma_electron_max_array[4], - "OP ", - ) - po.ovarre( - self.outfile, - "Hugill-Murakami Mq model", - "(nd_plasma_electron_max_array(6))", - physics_variables.nd_plasma_electron_max_array[5], - "OP ", - ) - po.ovarre( - self.outfile, - "Greenwald model", - "(nd_plasma_electron_max_array(7))", - physics_variables.nd_plasma_electron_max_array[6], - "OP ", - ) - po.ovarre( - self.outfile, - "ASDEX New", - "(nd_plasma_electron_max_array(8))", - physics_variables.nd_plasma_electron_max_array[7], - "OP ", - ) - po.ovarre( - self.outfile, - "Density limit from scaling (/m3)", - "(nd_plasma_electrons_max)", - physics_variables.nd_plasma_electrons_max, - "OP ", - ) - - po.oblnkl(self.outfile) - po.ostars(self.outfile, 110) - po.oblnkl(self.outfile) + self.density_limit.output_density_limit_information() po.oheadr(self.outfile, "Plasma Reactions :") @@ -10270,6 +10210,78 @@ def calculate_density_limit( i_density_limit - 1 ] + def output_density_limit_information(self): + """Output density limit information to file.""" + + po.osubhd(self.outfile, "Density Limit using different models :") + po.ovarre( + self.outfile, + "Old ASDEX model", + "(nd_plasma_electron_max_array(1))", + physics_variables.nd_plasma_electron_max_array[0], + "OP ", + ) + po.ovarre( + self.outfile, + "Borrass ITER model I", + "(nd_plasma_electron_max_array(2))", + physics_variables.nd_plasma_electron_max_array[1], + "OP ", + ) + po.ovarre( + self.outfile, + "Borrass ITER model II", + "(nd_plasma_electron_max_array(3))", + physics_variables.nd_plasma_electron_max_array[2], + "OP ", + ) + po.ovarre( + self.outfile, + "JET edge radiation model", + "(nd_plasma_electron_max_array(4))", + physics_variables.nd_plasma_electron_max_array[3], + "OP ", + ) + po.ovarre( + self.outfile, + "JET simplified model", + "(nd_plasma_electron_max_array(5))", + physics_variables.nd_plasma_electron_max_array[4], + "OP ", + ) + po.ovarre( + self.outfile, + "Hugill-Murakami Mq model", + "(nd_plasma_electron_max_array(6))", + physics_variables.nd_plasma_electron_max_array[5], + "OP ", + ) + po.ovarre( + self.outfile, + "Greenwald model", + "(nd_plasma_electron_max_array(7))", + physics_variables.nd_plasma_electron_max_array[6], + "OP ", + ) + po.ovarre( + self.outfile, + "ASDEX New", + "(nd_plasma_electron_max_array(8))", + physics_variables.nd_plasma_electron_max_array[7], + "OP ", + ) + po.ovarre( + self.outfile, + "Density limit from scaling (/m3)", + "(nd_plasma_electrons_max)", + physics_variables.nd_plasma_electrons_max, + "OP ", + ) + + po.oblnkl(self.outfile) + po.ostars(self.outfile, 110) + po.oblnkl(self.outfile) + class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" From fffa17a69faf59b28681314e359476ae51fd828a Mon Sep 17 00:00:00 2001 From: mn3981 Date: Sat, 21 Feb 2026 14:21:28 +0000 Subject: [PATCH 13/19] Post merge conflict rebase commit --- process/core/io/plot_proc.py | 3 --- process/main.py | 1 + tests/unit/test_physics.py | 2 +- tests/unit/test_stellarator.py | 7 ++++++- 4 files changed, 8 insertions(+), 5 deletions(-) diff --git a/process/core/io/plot_proc.py b/process/core/io/plot_proc.py index 42788f4c8..af8b0dd94 100644 --- a/process/core/io/plot_proc.py +++ b/process/core/io/plot_proc.py @@ -252,9 +252,6 @@ def plot_plasma( axis.plot(pg.rs[0], pg.zs[0], color="black") axis.plot(pg.rs[1], pg.zs[1], color="black") - axis.plot(pg.rs[0] + 0.005, pg.zs[0] + 0.005, color="red", linestyle="dashed") - axis.plot(pg.rs[1] + 0.005, pg.zs[1] + 0.005, color="red", linestyle="dashed") - # Set triang_95 to stop plotting plasma past boundary # Assume IPDG scaling triang_95 = triang / 1.5 diff --git a/process/main.py b/process/main.py index 9b504a21b..52f6b62b4 100644 --- a/process/main.py +++ b/process/main.py @@ -97,6 +97,7 @@ DetailedPhysics, Physics, PlasmaBeta, + PlasmaDensityLimit, PlasmaInductance, ) from process.models.physics.plasma_geometry import PlasmaGeom diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 3b5f7a5c9..f7c141787 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -24,8 +24,8 @@ DetailedPhysics, Physics, PlasmaBeta, - PlasmaInductance, PlasmaDensityLimit, + PlasmaInductance, calculate_current_coefficient_hastie, calculate_plasma_current_peng, calculate_poloidal_field, diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 1cd6fca15..9c641cd9e 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -29,7 +29,12 @@ LowerHybrid, NeutralBeam, ) -from process.models.physics.physics import Physics, PlasmaBeta, PlasmaInductance +from process.models.physics.physics import ( + Physics, + PlasmaBeta, + PlasmaDensityLimit, + PlasmaInductance, +) from process.models.physics.plasma_profiles import PlasmaProfile from process.models.power import Power from process.models.stellarator.build import st_build From 94465bf1582e1fc6a2f996ef87b6c03dec55a0df Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 26 Feb 2026 08:49:15 +0000 Subject: [PATCH 14/19] Requested changes, create new density limit file in physics models --- .../source/physics-models/plasma_density.md | 2 +- process/main.py | 2 +- process/models/physics/density_limit.py | 647 ++++++++++++++ process/models/physics/physics.py | 818 ------------------ tests/unit/test_physics.py | 2 +- tests/unit/test_stellarator.py | 2 +- 6 files changed, 651 insertions(+), 822 deletions(-) create mode 100644 process/models/physics/density_limit.py diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index 04f6ca812..0cdeae7ae 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -1,4 +1,4 @@ -# Density Limit +# Density Limit | `PlasmaDensityLimit` Several density limit models are available in PROCESS. These are calculated in routine `calculate_density_limit()`, which is called by `physics`. diff --git a/process/main.py b/process/main.py index 52f6b62b4..7388e750c 100644 --- a/process/main.py +++ b/process/main.py @@ -92,12 +92,12 @@ LowerHybrid, NeutralBeam, ) +from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, Physics, PlasmaBeta, - PlasmaDensityLimit, PlasmaInductance, ) from process.models.physics.plasma_geometry import PlasmaGeom diff --git a/process/models/physics/density_limit.py b/process/models/physics/density_limit.py new file mode 100644 index 000000000..d0aba41a8 --- /dev/null +++ b/process/models/physics/density_limit.py @@ -0,0 +1,647 @@ +import logging +from enum import IntEnum + +import numpy as np + +from process import constants +from process import process_output as po +from process.data_structure import ( + current_drive_variables, + divertor_variables, + physics_variables, +) +from process.exceptions import ProcessValueError + +logger = logging.getLogger(__name__) + + +class DensityLimitModel(IntEnum): + """Electron density model types""" + + ASDEX = 1 + BORRASS_ITER_I = 2 + BORRASS_ITER_II = 3 + JET_EDGE_RADIATION = 4 + JET_SIMPLE = 5 + HUGILL_MURAKAMI = 6 + GREENWALD = 7 + ASDEX_NEW = 8 + + +class PlasmaDensityLimit: + """Class to hold plasma density limit calculations for plasma processing.""" + + def __init__(self): + self.outfile = constants.NOUT + self.mfile = constants.MFILE + + def run(self): + nd_plasma_electron_max_array = np.empty((8,)) + + p_perp = ( + physics_variables.p_plasma_separatrix_mw / physics_variables.a_plasma_surface + ) + + # Old ASDEX density limit formula + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # Borrass density limit model for ITER (I) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # Borrass et al, ITER-TN-PH-9-6 (1989) + + nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # Borrass density limit model for ITER (II) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # This formula is (almost) identical to that in the original routine + # denlim (now deleted). + + nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + q95=physics_variables.q95, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # JET edge radiation density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # qcyl=qstar here, but literature is not clear. + + nd_plasma_electron_max_array[3] = ( + self.calculate_jet_edge_radiation_density_limit( + zeff=physics_variables.n_charge_plasma_effective_vol_avg, + p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, + prn1=divertor_variables.prn1, + qcyl=physics_variables.qstar, + ) + ) + + # JET simplified density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + p_plasma_separatrix_mw=physics_variables.p_plasma_separatrix_mw, + rmajor=physics_variables.rmajor, + prn1=divertor_variables.prn1, + ) + + # Hugill-Murakami M.q limit + # qcyl=qstar here, which is okay according to the literature + + nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + rmajor=physics_variables.rmajor, + qcyl=physics_variables.qstar, + ) + + # Greenwald limit + + nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( + c_plasma=physics_variables.plasma_current, rminor=physics_variables.rminor + ) + + nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( + p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, + c_plasma=physics_variables.plasma_current, + q95=physics_variables.q95, + prn1=divertor_variables.prn1, + ) + + physics_variables.nd_plasma_electron_max_array = nd_plasma_electron_max_array + + # Calculate beta_norm_max based on i_beta_norm_max + try: + model = DensityLimitModel(int(physics_variables.i_density_limit)) + physics_variables.nd_plasma_electrons_max = self.get_density_limit_value( + model + ) + except ValueError: + raise ProcessValueError( + "Illegal value of i_density_limit", + i_density_limit=physics_variables.i_density_limit, + ) from None + + def get_density_limit_value(self, model: DensityLimitModel) -> float: + """Get the density limit value (n_e_max) for the specified model.""" + model_map = { + DensityLimitModel.ASDEX: physics_variables.nd_plasma_electron_max_array[0], + DensityLimitModel.BORRASS_ITER_I: physics_variables.nd_plasma_electron_max_array[ + 1 + ], + DensityLimitModel.BORRASS_ITER_II: physics_variables.nd_plasma_electron_max_array[ + 2 + ], + DensityLimitModel.JET_EDGE_RADIATION: physics_variables.nd_plasma_electron_max_array[ + 3 + ], + DensityLimitModel.JET_SIMPLE: physics_variables.nd_plasma_electron_max_array[ + 4 + ], + DensityLimitModel.HUGILL_MURAKAMI: physics_variables.nd_plasma_electron_max_array[ + 5 + ], + DensityLimitModel.GREENWALD: physics_variables.nd_plasma_electron_max_array[ + 6 + ], + DensityLimitModel.ASDEX_NEW: physics_variables.nd_plasma_electron_max_array[ + 7 + ], + } + return model_map[model] + + @staticmethod + def calculate_asdex_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the ASDEX density limit. + + :param p_perp: Perpendicular power density (MW/m²). + :type p_perp: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param q95: Safety factor at 95% of the plasma poloidal flux. + :type q95: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The ASDEX density limit (m⁻³). + :rtype: float + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 1.54e20 + * p_perp**0.43 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.45 + ) / prn1 + + @staticmethod + def calculate_borrass_iter_i_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the Borrass ITER I density limit. + + :param p_perp: Perpendicular power density (MW/m²). + :type p_perp: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param q95: Safety factor at 95% of the plasma poloidal flux. + :type q95: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The Borrass ITER I density limit (m⁻³). + :rtype: float + + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 1.8e20 + * p_perp**0.53 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.22 + ) / prn1 + + @staticmethod + def calculate_borrass_iter_ii_density_limit( + p_perp: float, + b_plasma_toroidal_on_axis: float, + q95: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the Borrass ITER II density limit. + + :param p_perp: Perpendicular power density (MW/m²). + :type p_perp: float + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param q95: Safety factor at 95% of the plasma poloidal flux. + :type q95: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The Borrass ITER II density limit (m⁻³). + :rtype: float + + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 0.5e20 + * p_perp**0.57 + * b_plasma_toroidal_on_axis**0.31 + / (q95 * rmajor) ** 0.09 + ) / prn1 + + @staticmethod + def calculate_jet_edge_radiation_density_limit( + zeff: float, p_hcd_injected_total_mw: float, prn1: float, qcyl: float + ) -> float: + """ + Calculate the JET edge radiation density limit. + + :param zeff: Effective charge (Z_eff). + :type zeff: float + :param p_hcd_injected_total_mw: Power injected into the plasma (MW). + :type p_hcd_injected_total_mw: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :param qcyl: Equivalent cylindrical safety factor (qstar). + :type qcyl: float + :return: The JET edge radiation density limit (m⁻³). + :rtype: float + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + """ + + denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) + if denom <= 0.0: + return 0.0 + return (1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom)) / prn1 + + @staticmethod + def calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis: float, + p_plasma_separatrix_mw: float, + rmajor: float, + prn1: float, + ) -> float: + """ + Calculate the JET simple density limit. + + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param p_plasma_separatrix_mw: Power crossing the separatrix (MW). + :type p_plasma_separatrix_mw: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The JET simple density limit (m⁻³). + :rtype: float + + :references: + - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + + """ + return ( + 0.237e20 + * b_plasma_toroidal_on_axis + * np.sqrt(p_plasma_separatrix_mw) + / rmajor + ) / prn1 + + @staticmethod + def calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis: float, rmajor: float, qcyl: float + ) -> float: + """ + Calculate the Hugill-Murakami density limit. + + :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). + :type b_plasma_toroidal_on_axis: float + :param rmajor: Plasma major radius (m). + :type rmajor: float + :param qcyl: Equivalent cylindrical safety factor (qstar). + :type qcyl: float + :return: The Hugill-Murakami density limit (m⁻³). + :rtype: float + + :references: + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', + """ + + return 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) + + @staticmethod + def calculate_greenwald_density_limit(c_plasma: float, rminor: float) -> float: + """ + Calculate the Greenwald density limit (n_GW). + + :param c_plasma: Plasma current (A). + :type c_plasma: float + :param rminor: Plasma minor radius (m). + :type rminor: float + :return: The Greenwald density limit (m⁻³). + :rtype: float + + :notes: The Greenwald limit is typically applied to the line averaged electron density + + :references: + - M. Greenwald et al., “A new look at density limits in tokamaks,” + Nuclear Fusion, vol. 28, no. 12, pp. 2199-2207, Dec. 1988, + doi: https://doi.org/10.1088/0029-5515/28/12/009. + + - M. Greenwald, “Density limits in toroidal plasmas,” + Plasma Physics and Controlled Fusion, vol. 44, no. 8, pp. R27-R53, Jul. 2002, + doi: https://doi.org/10.1088/0741-3335/44/8/201. + """ + + return 1.0e14 * c_plasma / (np.pi * rminor**2) + + @staticmethod + def calculate_asdex_new_density_limit( + p_hcd_injected_total_mw: float, c_plasma: float, q95: float, prn1: float + ) -> float: + """ + Calculate the ASDEX Upgrade new density limit. + + :param p_hcd_injected_total_mw: Power injected into the plasma (MW). + :type p_hcd_injected_total_mw: float + :param plasma_current: Plasma current (A). + :type plasma_current: float + :param q95: Safety factor at 95% surface. + :type q95: float + :param prn1: Edge density / average plasma density. + :type prn1: float + :return: The ASDEX Upgrade new density limit (m⁻³). + :rtype: float + + :notes: This limit is for the separatrix density so wee scale by `prn1` to get it as a volume average + + :references: + + - J. W. Berkery et al., “Density limits as disruption forecasters for spherical tokamaks,” + Plasma Physics and Controlled Fusion, vol. 65, no. 9, pp. 095003-095003, Jul. 2023, + doi: https://doi.org/10.1088/1361-6587/ace476. + + - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” vol. 57, no. 1, pp. 014038-014038, Nov. 2014, + doi: https://doi.org/10.1088/0741-3335/57/1/014038. + """ + return ( + 1.0e20 + * 0.506 + * (p_hcd_injected_total_mw**0.396 * (c_plasma / 1.0e6) ** 0.265) + / (q95**0.323) + ) / prn1 + + def calculate_density_limit( + self, + b_plasma_toroidal_on_axis: float, + i_density_limit: int, + p_plasma_separatrix_mw: float, + p_hcd_injected_total_mw: float, + plasma_current: float, + prn1: float, + qcyl: float, + q95: float, + rmajor: float, + rminor: float, + a_plasma_surface: float, + zeff: float, + ) -> tuple[np.ndarray, float]: + """ + Calculate the density limit using various models. + + Args: + b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). + i_density_limit (int): Switch denoting which formula to enforce. + p_plasma_separatrix_mw (float): Power flowing to the edge plasma via charged particles (MW). + p_hcd_injected_total_mw (float): Power injected into the plasma (MW). + plasma_current (float): Plasma current (A). + prn1 (float): Edge density / average plasma density. + qcyl (float): Equivalent cylindrical safety factor (qstar). + q95 (float): Safety factor at 95% surface. + rmajor (float): Plasma major radius (m). + rminor (float): Plasma minor radius (m). + a_plasma_surface (float): Plasma surface area (m^2). + zeff (float): Plasma effective charge. + + Returns: + Tuple[np.ndarray, float]: A tuple containing: + - nd_plasma_electron_max_array (np.ndarray): Average plasma density limit using seven different models (m^-3). + - nd_plasma_electrons_max (float): Enforced average plasma density limit (m^-3). + + Raises: + ValueError: If i_density_limit is not between 1 and 7. + + Notes: + This routine calculates several different formulae for the density limit and enforces the one chosen by the user. + For i_density_limit = 1-5, 8, we scale the sepatrix density limit output by the ratio of the separatrix to volume averaged density + + References: + - AEA FUS 172: Physics Assessment for the European Reactor Study + + - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989 + + - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” + vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. ‌ + """ + + if i_density_limit < 1 or i_density_limit > 7: + raise ProcessValueError( + "Illegal value for i_density_limit", i_density_limit=i_density_limit + ) + + nd_plasma_electron_max_array = np.empty((8,)) + + # Power per unit area crossing the plasma edge + # (excludes radiation and neutrons) + + p_perp = p_plasma_separatrix_mw / a_plasma_surface + + # Old ASDEX density limit formula + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) + + # Borrass density limit model for ITER (I) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # Borrass et al, ITER-TN-PH-9-6 (1989) + + nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) + + # Borrass density limit model for ITER (II) + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # This formula is (almost) identical to that in the original routine + # denlim (now deleted). + + nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( + p_perp=p_perp, + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + q95=q95, + rmajor=rmajor, + prn1=prn1, + ) + + # JET edge radiation density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + # qcyl=qstar here, but literature is not clear. + + nd_plasma_electron_max_array[3] = ( + self.calculate_jet_edge_radiation_density_limit( + zeff=zeff, + p_hcd_injected_total_mw=p_hcd_injected_total_mw, + prn1=prn1, + qcyl=qcyl, + ) + ) + + # JET simplified density limit model + # This applies to the density at the plasma edge, so must be scaled + # to give the density limit applying to the average plasma density. + + nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, + p_plasma_separatrix_mw=p_plasma_separatrix_mw, + rmajor=rmajor, + prn1=prn1, + ) + + # Hugill-Murakami M.q limit + # qcyl=qstar here, which is okay according to the literature + + nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( + b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, rmajor=rmajor, qcyl=qcyl + ) + + # Greenwald limit + + nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( + c_plasma=plasma_current, rminor=rminor + ) + + nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( + p_hcd_injected_total_mw=p_hcd_injected_total_mw, + c_plasma=plasma_current, + q95=q95, + prn1=prn1, + ) + + # Enforce the chosen density limit + + return nd_plasma_electron_max_array, nd_plasma_electron_max_array[ + i_density_limit - 1 + ] + + def output_density_limit_information(self): + """Output density limit information to file.""" + + po.osubhd(self.outfile, "Density Limit using different models :") + po.ovarre( + self.outfile, + "Old ASDEX model", + "(nd_plasma_electron_max_array(1))", + physics_variables.nd_plasma_electron_max_array[0], + "OP ", + ) + po.ovarre( + self.outfile, + "Borrass ITER model I", + "(nd_plasma_electron_max_array(2))", + physics_variables.nd_plasma_electron_max_array[1], + "OP ", + ) + po.ovarre( + self.outfile, + "Borrass ITER model II", + "(nd_plasma_electron_max_array(3))", + physics_variables.nd_plasma_electron_max_array[2], + "OP ", + ) + po.ovarre( + self.outfile, + "JET edge radiation model", + "(nd_plasma_electron_max_array(4))", + physics_variables.nd_plasma_electron_max_array[3], + "OP ", + ) + po.ovarre( + self.outfile, + "JET simplified model", + "(nd_plasma_electron_max_array(5))", + physics_variables.nd_plasma_electron_max_array[4], + "OP ", + ) + po.ovarre( + self.outfile, + "Hugill-Murakami Mq model", + "(nd_plasma_electron_max_array(6))", + physics_variables.nd_plasma_electron_max_array[5], + "OP ", + ) + po.ovarre( + self.outfile, + "Greenwald model", + "(nd_plasma_electron_max_array(7))", + physics_variables.nd_plasma_electron_max_array[6], + "OP ", + ) + po.ovarre( + self.outfile, + "ASDEX New", + "(nd_plasma_electron_max_array(8))", + physics_variables.nd_plasma_electron_max_array[7], + "OP ", + ) + po.ovarre( + self.outfile, + "Density limit from scaling (/m3)", + "(nd_plasma_electrons_max)", + physics_variables.nd_plasma_electrons_max, + "OP ", + ) + + po.oblnkl(self.outfile) + po.ostars(self.outfile, 110) + po.oblnkl(self.outfile) diff --git a/process/models/physics/physics.py b/process/models/physics/physics.py index db7ff7c20..02e925af3 100644 --- a/process/models/physics/physics.py +++ b/process/models/physics/physics.py @@ -2688,192 +2688,6 @@ def calculate_current_profile_index_wesson(qstar: float, q0: float) -> float: """ return qstar / q0 - 1.0 - @staticmethod - def calculate_density_limit( - b_plasma_toroidal_on_axis: float, - i_density_limit: int, - p_plasma_separatrix_mw: float, - p_hcd_injected_total_mw: float, - plasma_current: float, - prn1: float, - qcyl: float, - q95: float, - rmajor: float, - rminor: float, - a_plasma_surface: float, - zeff: float, - ) -> tuple[np.ndarray, float]: - """Calculate the density limit using various models. - - Parameters - ---------- - b_plasma_toroidal_on_axis : float - Toroidal field on axis (T). - i_density_limit : int - Switch denoting which formula to enforce. - p_plasma_separatrix_mw : float - Power flowing to the edge plasma via charged particles (MW). - p_hcd_injected_total_mw : float - Power injected into the plasma (MW). - plasma_current : float - Plasma current (A). - prn1 : float - Edge density / average plasma density. - qcyl : float - Equivalent cylindrical safety factor (qstar). - q95 : float - Safety factor at 95% surface. - rmajor : float - Plasma major radius (m). - rminor : float - Plasma minor radius (m). - a_plasma_surface : float - Plasma surface area (m^2). - zeff : float - Plasma effective charge. - - Returns - ------- - Tuple[np.ndarray, float] - A tuple containing: - - nd_plasma_electron_max_array (np.ndarray): Average plasma density limit using seven different models (m^-3). - - nd_plasma_electrons_max (float): Enforced average plasma density limit (m^-3). - - Raises - ------ - ValueError - If i_density_limit is not between 1 and 7. - Notes - ----- - ValueError - If i_density_limit is not between 1 and 7. - Notes - ----- - This routine calculates several different formulae for the density limit and enforces the one chosen by the user. - ValueError - If i_density_limit is not between 1 and 7. - Notes - ----- - This routine calculates several different formulae for the density limit and enforces the one chosen by the user. - For i_density_limit = 1-5, 8, we scale the sepatrix density limit output by the ratio of the separatrix to volume averaged density - - References - ---------- - - AEA FUS 172 - Physics Assessment for the European Reactor Study - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines - 1989 - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines - 1989 - - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” - vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. ‌ - - """ - - if i_density_limit < 1 or i_density_limit > 7: - raise ProcessValueError( - "Illegal value for i_density_limit", i_density_limit=i_density_limit - ) - - nd_plasma_electron_max_array = np.empty((8,)) - - # Power per unit area crossing the plasma edge - # (excludes radiation and neutrons) - - p_perp = p_plasma_separatrix_mw / a_plasma_surface - - # Old ASDEX density limit formula - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[0] = ( - 1.54e20 - * p_perp**0.43 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.45 - ) / prn1 - - # Borrass density limit model for ITER (I) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # Borrass et al, ITER-TN-PH-9-6 (1989) - - nd_plasma_electron_max_array[1] = ( - 1.8e20 - * p_perp**0.53 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.22 - ) / prn1 - - # Borrass density limit model for ITER (II) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # This formula is (almost) identical to that in the original routine - # denlim (now deleted). - - nd_plasma_electron_max_array[2] = ( - 0.5e20 - * p_perp**0.57 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.09 - ) / prn1 - - # JET edge radiation density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # qcyl=qstar here, but literature is not clear. - - denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) - if denom <= 0.0: - if i_density_limit == 4: - logger.error( - f"qcyl < 4/3; nd_plasma_electron_max_array(4) set to zero; model 5 will be enforced instead. {denom=} {qcyl=}" - ) - i_density_limit = 5 - - nd_plasma_electron_max_array[3] = 0.0 - else: - nd_plasma_electron_max_array[3] = ( - 1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom) - ) / prn1 - - # JET simplified density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[4] = ( - 0.237e20 - * b_plasma_toroidal_on_axis - * np.sqrt(p_plasma_separatrix_mw) - / rmajor - ) / prn1 - - # Hugill-Murakami M.q limit - # qcyl=qstar here, which is okay according to the literature - - nd_plasma_electron_max_array[5] = ( - 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) - ) - - # Greenwald limit - - nd_plasma_electron_max_array[6] = ( - 1.0e14 * plasma_current / (np.pi * rminor * rminor) - ) - - nd_plasma_electron_max_array[7] = ( - 1.0e20 - * 0.506 - * (p_hcd_injected_total_mw**0.396 * (plasma_current / 1.0e6) ** 0.265) - / (q95**0.323) - ) / prn1 - - # Enforce the chosen density limit - - return nd_plasma_electron_max_array, nd_plasma_electron_max_array[ - i_density_limit - 1 - ] - @staticmethod def plasma_composition(): """Calculates various plasma component fractional makeups. @@ -9651,638 +9465,6 @@ def output_volt_second_information(self): po.oblnkl(self.outfile) -class DensityLimitModel(IntEnum): - """Electron density model types""" - - ASDEX = 1 - BORRASS_ITER_I = 2 - BORRASS_ITER_II = 3 - JET_EDGE_RADIATION = 4 - JET_SIMPLE = 5 - HUGILL_MURAKAMI = 6 - GREENWALD = 7 - ASDEX_NEW = 8 - - -class PlasmaDensityLimit: - """Class to hold plasma density limit calculations for plasma processing.""" - - def __init__(self): - self.outfile = constants.NOUT - self.mfile = constants.MFILE - - def run(self): - nd_plasma_electron_max_array = np.empty((8,)) - - p_perp = ( - physics_variables.p_plasma_separatrix_mw / physics_variables.a_plasma_surface - ) - - # Old ASDEX density limit formula - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - q95=physics_variables.q95, - rmajor=physics_variables.rmajor, - prn1=divertor_variables.prn1, - ) - - # Borrass density limit model for ITER (I) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # Borrass et al, ITER-TN-PH-9-6 (1989) - - nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - q95=physics_variables.q95, - rmajor=physics_variables.rmajor, - prn1=divertor_variables.prn1, - ) - - # Borrass density limit model for ITER (II) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # This formula is (almost) identical to that in the original routine - # denlim (now deleted). - - nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - q95=physics_variables.q95, - rmajor=physics_variables.rmajor, - prn1=divertor_variables.prn1, - ) - - # JET edge radiation density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # qcyl=qstar here, but literature is not clear. - - nd_plasma_electron_max_array[3] = ( - self.calculate_jet_edge_radiation_density_limit( - zeff=physics_variables.n_charge_plasma_effective_vol_avg, - p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, - prn1=divertor_variables.prn1, - qcyl=physics_variables.qstar, - ) - ) - - # JET simplified density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - p_plasma_separatrix_mw=physics_variables.p_plasma_separatrix_mw, - rmajor=physics_variables.rmajor, - prn1=divertor_variables.prn1, - ) - - # Hugill-Murakami M.q limit - # qcyl=qstar here, which is okay according to the literature - - nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - rmajor=physics_variables.rmajor, - qcyl=physics_variables.qstar, - ) - - # Greenwald limit - - nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( - c_plasma=physics_variables.plasma_current, rminor=physics_variables.rminor - ) - - nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( - p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, - c_plasma=physics_variables.plasma_current, - q95=physics_variables.q95, - prn1=divertor_variables.prn1, - ) - - physics_variables.nd_plasma_electron_max_array = nd_plasma_electron_max_array - - # Calculate beta_norm_max based on i_beta_norm_max - try: - model = DensityLimitModel(int(physics_variables.i_density_limit)) - physics_variables.nd_plasma_electrons_max = self.get_density_limit_value( - model - ) - except ValueError: - raise ProcessValueError( - "Illegal value of i_density_limit", - i_density_limit=physics_variables.i_density_limit, - ) from None - - def get_density_limit_value(self, model: DensityLimitModel) -> float: - """Get the density limit value (n_e_max) for the specified model.""" - model_map = { - DensityLimitModel.ASDEX: physics_variables.nd_plasma_electron_max_array[0], - DensityLimitModel.BORRASS_ITER_I: physics_variables.nd_plasma_electron_max_array[ - 1 - ], - DensityLimitModel.BORRASS_ITER_II: physics_variables.nd_plasma_electron_max_array[ - 2 - ], - DensityLimitModel.JET_EDGE_RADIATION: physics_variables.nd_plasma_electron_max_array[ - 3 - ], - DensityLimitModel.JET_SIMPLE: physics_variables.nd_plasma_electron_max_array[ - 4 - ], - DensityLimitModel.HUGILL_MURAKAMI: physics_variables.nd_plasma_electron_max_array[ - 5 - ], - DensityLimitModel.GREENWALD: physics_variables.nd_plasma_electron_max_array[ - 6 - ], - DensityLimitModel.ASDEX_NEW: physics_variables.nd_plasma_electron_max_array[ - 7 - ], - } - return model_map[model] - - @staticmethod - def calculate_asdex_density_limit( - p_perp: float, - b_plasma_toroidal_on_axis: float, - q95: float, - rmajor: float, - prn1: float, - ) -> float: - """ - Calculate the ASDEX density limit. - - :param p_perp: Perpendicular power density (MW/m²). - :type p_perp: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param q95: Safety factor at 95% of the plasma poloidal flux. - :type q95: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The ASDEX density limit (m⁻³). - :rtype: float - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - return ( - 1.54e20 - * p_perp**0.43 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.45 - ) / prn1 - - @staticmethod - def calculate_borrass_iter_i_density_limit( - p_perp: float, - b_plasma_toroidal_on_axis: float, - q95: float, - rmajor: float, - prn1: float, - ) -> float: - """ - Calculate the Borrass ITER I density limit. - - :param p_perp: Perpendicular power density (MW/m²). - :type p_perp: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param q95: Safety factor at 95% of the plasma poloidal flux. - :type q95: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The Borrass ITER I density limit (m⁻³). - :rtype: float - - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - return ( - 1.8e20 - * p_perp**0.53 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.22 - ) / prn1 - - @staticmethod - def calculate_borrass_iter_ii_density_limit( - p_perp: float, - b_plasma_toroidal_on_axis: float, - q95: float, - rmajor: float, - prn1: float, - ) -> float: - """ - Calculate the Borrass ITER II density limit. - - :param p_perp: Perpendicular power density (MW/m²). - :type p_perp: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param q95: Safety factor at 95% of the plasma poloidal flux. - :type q95: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The Borrass ITER II density limit (m⁻³). - :rtype: float - - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - return ( - 0.5e20 - * p_perp**0.57 - * b_plasma_toroidal_on_axis**0.31 - / (q95 * rmajor) ** 0.09 - ) / prn1 - - @staticmethod - def calculate_jet_edge_radiation_density_limit( - zeff: float, p_hcd_injected_total_mw: float, prn1: float, qcyl: float - ) -> float: - """ - Calculate the JET edge radiation density limit. - - :param zeff: Effective charge (Z_eff). - :type zeff: float - :param p_hcd_injected_total_mw: Power injected into the plasma (MW). - :type p_hcd_injected_total_mw: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :param qcyl: Equivalent cylindrical safety factor (qstar). - :type qcyl: float - :return: The JET edge radiation density limit (m⁻³). - :rtype: float - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - """ - - denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) - if denom <= 0.0: - return 0.0 - return (1.0e20 * np.sqrt(p_hcd_injected_total_mw / denom)) / prn1 - - @staticmethod - def calculate_jet_simple_density_limit( - b_plasma_toroidal_on_axis: float, - p_plasma_separatrix_mw: float, - rmajor: float, - prn1: float, - ) -> float: - """ - Calculate the JET simple density limit. - - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param p_plasma_separatrix_mw: Power crossing the separatrix (MW). - :type p_plasma_separatrix_mw: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The JET simple density limit (m⁻³). - :rtype: float - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - - """ - return ( - 0.237e20 - * b_plasma_toroidal_on_axis - * np.sqrt(p_plasma_separatrix_mw) - / rmajor - ) / prn1 - - @staticmethod - def calculate_hugill_murakami_density_limit( - b_plasma_toroidal_on_axis: float, rmajor: float, qcyl: float - ) -> float: - """ - Calculate the Hugill-Murakami density limit. - - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param qcyl: Equivalent cylindrical safety factor (qstar). - :type qcyl: float - :return: The Hugill-Murakami density limit (m⁻³). - :rtype: float - - :references: - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', - """ - - return 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) - - @staticmethod - def calculate_greenwald_density_limit(c_plasma: float, rminor: float) -> float: - """ - Calculate the Greenwald density limit (n_GW). - - :param c_plasma: Plasma current (A). - :type c_plasma: float - :param rminor: Plasma minor radius (m). - :type rminor: float - :return: The Greenwald density limit (m⁻³). - :rtype: float - - :notes: The Greenwald limit is typically applied to the line averaged electron density - - :references: - - M. Greenwald et al., “A new look at density limits in tokamaks,” - Nuclear Fusion, vol. 28, no. 12, pp. 2199-2207, Dec. 1988, - doi: https://doi.org/10.1088/0029-5515/28/12/009. - - - M. Greenwald, “Density limits in toroidal plasmas,” - Plasma Physics and Controlled Fusion, vol. 44, no. 8, pp. R27-R53, Jul. 2002, - doi: https://doi.org/10.1088/0741-3335/44/8/201. - """ - - return 1.0e14 * c_plasma / (np.pi * rminor**2) - - @staticmethod - def calculate_asdex_new_density_limit( - p_hcd_injected_total_mw: float, c_plasma: float, q95: float, prn1: float - ) -> float: - """ - Calculate the ASDEX Upgrade new density limit. - - :param p_hcd_injected_total_mw: Power injected into the plasma (MW). - :type p_hcd_injected_total_mw: float - :param plasma_current: Plasma current (A). - :type plasma_current: float - :param q95: Safety factor at 95% surface. - :type q95: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The ASDEX Upgrade new density limit (m⁻³). - :rtype: float - - :notes: This limit is for the separatrix density so wee scale by `prn1` to get it as a volume average - - :references: - - - J. W. Berkery et al., “Density limits as disruption forecasters for spherical tokamaks,” - Plasma Physics and Controlled Fusion, vol. 65, no. 9, pp. 095003-095003, Jul. 2023, - doi: https://doi.org/10.1088/1361-6587/ace476. - - - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” vol. 57, no. 1, pp. 014038-014038, Nov. 2014, - doi: https://doi.org/10.1088/0741-3335/57/1/014038. - """ - return ( - 1.0e20 - * 0.506 - * (p_hcd_injected_total_mw**0.396 * (c_plasma / 1.0e6) ** 0.265) - / (q95**0.323) - ) / prn1 - - def calculate_density_limit( - self, - b_plasma_toroidal_on_axis: float, - i_density_limit: int, - p_plasma_separatrix_mw: float, - p_hcd_injected_total_mw: float, - plasma_current: float, - prn1: float, - qcyl: float, - q95: float, - rmajor: float, - rminor: float, - a_plasma_surface: float, - zeff: float, - ) -> tuple[np.ndarray, float]: - """ - Calculate the density limit using various models. - - Args: - b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). - i_density_limit (int): Switch denoting which formula to enforce. - p_plasma_separatrix_mw (float): Power flowing to the edge plasma via charged particles (MW). - p_hcd_injected_total_mw (float): Power injected into the plasma (MW). - plasma_current (float): Plasma current (A). - prn1 (float): Edge density / average plasma density. - qcyl (float): Equivalent cylindrical safety factor (qstar). - q95 (float): Safety factor at 95% surface. - rmajor (float): Plasma major radius (m). - rminor (float): Plasma minor radius (m). - a_plasma_surface (float): Plasma surface area (m^2). - zeff (float): Plasma effective charge. - - Returns: - Tuple[np.ndarray, float]: A tuple containing: - - nd_plasma_electron_max_array (np.ndarray): Average plasma density limit using seven different models (m^-3). - - nd_plasma_electrons_max (float): Enforced average plasma density limit (m^-3). - - Raises: - ValueError: If i_density_limit is not between 1 and 7. - - Notes: - This routine calculates several different formulae for the density limit and enforces the one chosen by the user. - For i_density_limit = 1-5, 8, we scale the sepatrix density limit output by the ratio of the separatrix to volume averaged density - - References: - - AEA FUS 172: Physics Assessment for the European Reactor Study - - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989 - - - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” - vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. ‌ - """ - - if i_density_limit < 1 or i_density_limit > 7: - raise ProcessValueError( - "Illegal value for i_density_limit", i_density_limit=i_density_limit - ) - - nd_plasma_electron_max_array = np.empty((8,)) - - # Power per unit area crossing the plasma edge - # (excludes radiation and neutrons) - - p_perp = p_plasma_separatrix_mw / a_plasma_surface - - # Old ASDEX density limit formula - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, - q95=q95, - rmajor=rmajor, - prn1=prn1, - ) - - # Borrass density limit model for ITER (I) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # Borrass et al, ITER-TN-PH-9-6 (1989) - - nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, - q95=q95, - rmajor=rmajor, - prn1=prn1, - ) - - # Borrass density limit model for ITER (II) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # This formula is (almost) identical to that in the original routine - # denlim (now deleted). - - nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, - q95=q95, - rmajor=rmajor, - prn1=prn1, - ) - - # JET edge radiation density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # qcyl=qstar here, but literature is not clear. - - nd_plasma_electron_max_array[3] = ( - self.calculate_jet_edge_radiation_density_limit( - zeff=zeff, - p_hcd_injected_total_mw=p_hcd_injected_total_mw, - prn1=prn1, - qcyl=qcyl, - ) - ) - - # JET simplified density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( - b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, - p_plasma_separatrix_mw=p_plasma_separatrix_mw, - rmajor=rmajor, - prn1=prn1, - ) - - # Hugill-Murakami M.q limit - # qcyl=qstar here, which is okay according to the literature - - nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( - b_plasma_toroidal_on_axis=b_plasma_toroidal_on_axis, rmajor=rmajor, qcyl=qcyl - ) - - # Greenwald limit - - nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( - c_plasma=plasma_current, rminor=rminor - ) - - nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( - p_hcd_injected_total_mw=p_hcd_injected_total_mw, - c_plasma=plasma_current, - q95=q95, - prn1=prn1, - ) - - # Enforce the chosen density limit - - return nd_plasma_electron_max_array, nd_plasma_electron_max_array[ - i_density_limit - 1 - ] - - def output_density_limit_information(self): - """Output density limit information to file.""" - - po.osubhd(self.outfile, "Density Limit using different models :") - po.ovarre( - self.outfile, - "Old ASDEX model", - "(nd_plasma_electron_max_array(1))", - physics_variables.nd_plasma_electron_max_array[0], - "OP ", - ) - po.ovarre( - self.outfile, - "Borrass ITER model I", - "(nd_plasma_electron_max_array(2))", - physics_variables.nd_plasma_electron_max_array[1], - "OP ", - ) - po.ovarre( - self.outfile, - "Borrass ITER model II", - "(nd_plasma_electron_max_array(3))", - physics_variables.nd_plasma_electron_max_array[2], - "OP ", - ) - po.ovarre( - self.outfile, - "JET edge radiation model", - "(nd_plasma_electron_max_array(4))", - physics_variables.nd_plasma_electron_max_array[3], - "OP ", - ) - po.ovarre( - self.outfile, - "JET simplified model", - "(nd_plasma_electron_max_array(5))", - physics_variables.nd_plasma_electron_max_array[4], - "OP ", - ) - po.ovarre( - self.outfile, - "Hugill-Murakami Mq model", - "(nd_plasma_electron_max_array(6))", - physics_variables.nd_plasma_electron_max_array[5], - "OP ", - ) - po.ovarre( - self.outfile, - "Greenwald model", - "(nd_plasma_electron_max_array(7))", - physics_variables.nd_plasma_electron_max_array[6], - "OP ", - ) - po.ovarre( - self.outfile, - "ASDEX New", - "(nd_plasma_electron_max_array(8))", - physics_variables.nd_plasma_electron_max_array[7], - "OP ", - ) - po.ovarre( - self.outfile, - "Density limit from scaling (/m3)", - "(nd_plasma_electrons_max)", - physics_variables.nd_plasma_electrons_max, - "OP ", - ) - - po.oblnkl(self.outfile) - po.ostars(self.outfile, 110) - po.oblnkl(self.outfile) - - class DetailedPhysics: """Class to hold detailed physics models for plasma processing.""" diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index f7c141787..9a5d91097 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -19,12 +19,12 @@ LowerHybrid, NeutralBeam, ) +from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.impurity_radiation import initialise_imprad from process.models.physics.physics import ( DetailedPhysics, Physics, PlasmaBeta, - PlasmaDensityLimit, PlasmaInductance, calculate_current_coefficient_hastie, calculate_plasma_current_peng, diff --git a/tests/unit/test_stellarator.py b/tests/unit/test_stellarator.py index 9c641cd9e..3809b0fa8 100644 --- a/tests/unit/test_stellarator.py +++ b/tests/unit/test_stellarator.py @@ -29,10 +29,10 @@ LowerHybrid, NeutralBeam, ) +from process.models.physics.density_limit import PlasmaDensityLimit from process.models.physics.physics import ( Physics, PlasmaBeta, - PlasmaDensityLimit, PlasmaInductance, ) from process.models.physics.plasma_profiles import PlasmaProfile From cc5a99c116ebbd1add72bdc82c08d1a9db1d336f Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 26 Feb 2026 13:15:20 +0000 Subject: [PATCH 15/19] Fix documentation: update routine name for density limit calculation --- documentation/source/physics-models/plasma_density.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/documentation/source/physics-models/plasma_density.md b/documentation/source/physics-models/plasma_density.md index 0cdeae7ae..5b07350b9 100644 --- a/documentation/source/physics-models/plasma_density.md +++ b/documentation/source/physics-models/plasma_density.md @@ -1,7 +1,7 @@ # Density Limit | `PlasmaDensityLimit` Several density limit models are available in PROCESS. These are -calculated in routine `calculate_density_limit()`, which is called by `physics`. +calculated in the run routine `density_limit.run()`, which is called by `physics()`. This constraint can be activated by stating `icc = 5` in the input file. From a53c29d3a3b64fffb46d389fc1ef7878880ae58b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 26 Feb 2026 13:22:30 +0000 Subject: [PATCH 16/19] Change to numpy docstring style --- process/models/physics/density_limit.py | 396 ++++++++++++++---------- 1 file changed, 236 insertions(+), 160 deletions(-) diff --git a/process/models/physics/density_limit.py b/process/models/physics/density_limit.py index d0aba41a8..a3b0541aa 100644 --- a/process/models/physics/density_limit.py +++ b/process/models/physics/density_limit.py @@ -142,7 +142,19 @@ def run(self): ) from None def get_density_limit_value(self, model: DensityLimitModel) -> float: - """Get the density limit value (n_e_max) for the specified model.""" + """ + Get the density limit value (n_e_max) for the specified model. + + Parameters + ---------- + model : DensityLimitModel + The density limit model type. + + Returns + ------- + float + The density limit value (m⁻³). + """ model_map = { DensityLimitModel.ASDEX: physics_variables.nd_plasma_electron_max_array[0], DensityLimitModel.BORRASS_ITER_I: physics_variables.nd_plasma_electron_max_array[ @@ -180,22 +192,27 @@ def calculate_asdex_density_limit( """ Calculate the ASDEX density limit. - :param p_perp: Perpendicular power density (MW/m²). - :type p_perp: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param q95: Safety factor at 95% of the plasma poloidal flux. - :type q95: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The ASDEX density limit (m⁻³). - :rtype: float - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - + Parameters + ---------- + p_perp : float + Perpendicular power density (MW/m²). + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + q95 : float + Safety factor at 95% of the plasma poloidal flux. + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The ASDEX density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 """ return ( 1.54e20 @@ -215,23 +232,27 @@ def calculate_borrass_iter_i_density_limit( """ Calculate the Borrass ITER I density limit. - :param p_perp: Perpendicular power density (MW/m²). - :type p_perp: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param q95: Safety factor at 95% of the plasma poloidal flux. - :type q95: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The Borrass ITER I density limit (m⁻³). - :rtype: float - - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - + Parameters + ---------- + p_perp : float + Perpendicular power density (MW/m²). + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + q95 : float + Safety factor at 95% of the plasma poloidal flux. + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The Borrass ITER I density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 """ return ( 1.8e20 @@ -251,23 +272,27 @@ def calculate_borrass_iter_ii_density_limit( """ Calculate the Borrass ITER II density limit. - :param p_perp: Perpendicular power density (MW/m²). - :type p_perp: float - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param q95: Safety factor at 95% of the plasma poloidal flux. - :type q95: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The Borrass ITER II density limit (m⁻³). - :rtype: float - - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - + Parameters + ---------- + p_perp : float + Perpendicular power density (MW/m²). + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + q95 : float + Safety factor at 95% of the plasma poloidal flux. + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The Borrass ITER II density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 """ return ( 0.5e20 @@ -283,19 +308,25 @@ def calculate_jet_edge_radiation_density_limit( """ Calculate the JET edge radiation density limit. - :param zeff: Effective charge (Z_eff). - :type zeff: float - :param p_hcd_injected_total_mw: Power injected into the plasma (MW). - :type p_hcd_injected_total_mw: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :param qcyl: Equivalent cylindrical safety factor (qstar). - :type qcyl: float - :return: The JET edge radiation density limit (m⁻³). - :rtype: float - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 + Parameters + ---------- + zeff : float + Effective charge (Z_eff). + p_hcd_injected_total_mw : float + Power injected into the plasma (MW). + prn1 : float + Edge density / average plasma density. + qcyl : float + Equivalent cylindrical safety factor (qstar). + + Returns + ------- + float + The JET edge radiation density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 """ denom = (zeff - 1.0) * (1.0 - 4.0 / (3.0 * qcyl)) @@ -313,20 +344,25 @@ def calculate_jet_simple_density_limit( """ Calculate the JET simple density limit. - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param p_plasma_separatrix_mw: Power crossing the separatrix (MW). - :type p_plasma_separatrix_mw: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The JET simple density limit (m⁻³). - :rtype: float - - :references: - - T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 - + Parameters + ---------- + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + p_plasma_separatrix_mw : float + Power crossing the separatrix (MW). + rmajor : float + Plasma major radius (m). + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The JET simple density limit (m⁻³). + + References + ---------- + T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992 """ return ( 0.237e20 @@ -342,17 +378,23 @@ def calculate_hugill_murakami_density_limit( """ Calculate the Hugill-Murakami density limit. - :param b_plasma_toroidal_on_axis: Toroidal field on axis (T). - :type b_plasma_toroidal_on_axis: float - :param rmajor: Plasma major radius (m). - :type rmajor: float - :param qcyl: Equivalent cylindrical safety factor (qstar). - :type qcyl: float - :return: The Hugill-Murakami density limit (m⁻³). - :rtype: float - - :references: - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989', + Parameters + ---------- + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + rmajor : float + Plasma major radius (m). + qcyl : float + Equivalent cylindrical safety factor (qstar). + + Returns + ------- + float + The Hugill-Murakami density limit (m⁻³). + + References + ---------- + N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' """ return 3.0e20 * b_plasma_toroidal_on_axis / (rmajor * qcyl) @@ -362,23 +404,31 @@ def calculate_greenwald_density_limit(c_plasma: float, rminor: float) -> float: """ Calculate the Greenwald density limit (n_GW). - :param c_plasma: Plasma current (A). - :type c_plasma: float - :param rminor: Plasma minor radius (m). - :type rminor: float - :return: The Greenwald density limit (m⁻³). - :rtype: float - - :notes: The Greenwald limit is typically applied to the line averaged electron density - - :references: - - M. Greenwald et al., “A new look at density limits in tokamaks,” - Nuclear Fusion, vol. 28, no. 12, pp. 2199-2207, Dec. 1988, - doi: https://doi.org/10.1088/0029-5515/28/12/009. - - - M. Greenwald, “Density limits in toroidal plasmas,” - Plasma Physics and Controlled Fusion, vol. 44, no. 8, pp. R27-R53, Jul. 2002, - doi: https://doi.org/10.1088/0741-3335/44/8/201. + Parameters + ---------- + c_plasma : float + Plasma current (A). + rminor : float + Plasma minor radius (m). + + Returns + ------- + float + The Greenwald density limit (m⁻³). + + Notes + ----- + The Greenwald limit is typically applied to the line averaged electron density. + + References + ---------- + M. Greenwald et al., "A new look at density limits in tokamaks," + Nuclear Fusion, vol. 28, no. 12, pp. 2199-2207, Dec. 1988, + doi: https://doi.org/10.1088/0029-5515/28/12/009. + + M. Greenwald, "Density limits in toroidal plasmas," + Plasma Physics and Controlled Fusion, vol. 44, no. 8, pp. R27-R53, Jul. 2002, + doi: https://doi.org/10.1088/0741-3335/44/8/201. """ return 1.0e14 * c_plasma / (np.pi * rminor**2) @@ -390,27 +440,34 @@ def calculate_asdex_new_density_limit( """ Calculate the ASDEX Upgrade new density limit. - :param p_hcd_injected_total_mw: Power injected into the plasma (MW). - :type p_hcd_injected_total_mw: float - :param plasma_current: Plasma current (A). - :type plasma_current: float - :param q95: Safety factor at 95% surface. - :type q95: float - :param prn1: Edge density / average plasma density. - :type prn1: float - :return: The ASDEX Upgrade new density limit (m⁻³). - :rtype: float - - :notes: This limit is for the separatrix density so wee scale by `prn1` to get it as a volume average - - :references: - - - J. W. Berkery et al., “Density limits as disruption forecasters for spherical tokamaks,” - Plasma Physics and Controlled Fusion, vol. 65, no. 9, pp. 095003-095003, Jul. 2023, - doi: https://doi.org/10.1088/1361-6587/ace476. - - - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” vol. 57, no. 1, pp. 014038-014038, Nov. 2014, - doi: https://doi.org/10.1088/0741-3335/57/1/014038. + Parameters + ---------- + p_hcd_injected_total_mw : float + Power injected into the plasma (MW). + c_plasma : float + Plasma current (A). + q95 : float + Safety factor at 95% surface. + prn1 : float + Edge density / average plasma density. + + Returns + ------- + float + The ASDEX Upgrade new density limit (m⁻³). + + Notes + ----- + This limit is for the separatrix density so we scale by `prn1` to get it as a volume average. + + References + ---------- + J. W. Berkery et al., "Density limits as disruption forecasters for spherical tokamaks," + Plasma Physics and Controlled Fusion, vol. 65, no. 9, pp. 095003-095003, Jul. 2023, + doi: https://doi.org/10.1088/1361-6587/ace476. + + M. Bernert et al., "The H-mode density limit in the full tungsten ASDEX Upgrade tokamak," vol. 57, no. 1, pp. 014038-014038, Nov. 2014, + doi: https://doi.org/10.1088/0741-3335/57/1/014038. """ return ( 1.0e20 @@ -437,39 +494,58 @@ def calculate_density_limit( """ Calculate the density limit using various models. - Args: - b_plasma_toroidal_on_axis (float): Toroidal field on axis (T). - i_density_limit (int): Switch denoting which formula to enforce. - p_plasma_separatrix_mw (float): Power flowing to the edge plasma via charged particles (MW). - p_hcd_injected_total_mw (float): Power injected into the plasma (MW). - plasma_current (float): Plasma current (A). - prn1 (float): Edge density / average plasma density. - qcyl (float): Equivalent cylindrical safety factor (qstar). - q95 (float): Safety factor at 95% surface. - rmajor (float): Plasma major radius (m). - rminor (float): Plasma minor radius (m). - a_plasma_surface (float): Plasma surface area (m^2). - zeff (float): Plasma effective charge. - - Returns: - Tuple[np.ndarray, float]: A tuple containing: - - nd_plasma_electron_max_array (np.ndarray): Average plasma density limit using seven different models (m^-3). - - nd_plasma_electrons_max (float): Enforced average plasma density limit (m^-3). - - Raises: - ValueError: If i_density_limit is not between 1 and 7. - - Notes: - This routine calculates several different formulae for the density limit and enforces the one chosen by the user. - For i_density_limit = 1-5, 8, we scale the sepatrix density limit output by the ratio of the separatrix to volume averaged density - - References: - - AEA FUS 172: Physics Assessment for the European Reactor Study - - - N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989 - - - M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” - vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. ‌ + Parameters + ---------- + b_plasma_toroidal_on_axis : float + Toroidal field on axis (T). + i_density_limit : int + Switch denoting which formula to enforce (1-7). + p_plasma_separatrix_mw : float + Power flowing to the edge plasma via charged particles (MW). + p_hcd_injected_total_mw : float + Power injected into the plasma (MW). + plasma_current : float + Plasma current (A). + prn1 : float + Edge density / average plasma density. + qcyl : float + Equivalent cylindrical safety factor (qstar). + q95 : float + Safety factor at 95% surface. + rmajor : float + Plasma major radius (m). + rminor : float + Plasma minor radius (m). + a_plasma_surface : float + Plasma surface area (m²). + zeff : float + Plasma effective charge. + + Returns + ------- + tuple[np.ndarray, float] + A tuple containing: + - nd_plasma_electron_max_array : Average plasma density limit using eight different models (m⁻³). + - nd_plasma_electrons_max : Enforced average plasma density limit (m⁻³). + + Raises + ------ + ValueError + If i_density_limit is not between 1 and 7. + + Notes + ----- + This routine calculates several different formulae for the density limit and enforces the one chosen by the user. + For i_density_limit = 1-5, 8, we scale the separatrix density limit output by the ratio of the separatrix to volume averaged density. + + References + ---------- + AEA FUS 172: Physics Assessment for the European Reactor Study + + N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989' + + M. Bernert et al., "The H-mode density limit in the full tungsten ASDEX Upgrade tokamak," + vol. 57, no. 1, pp. 014038-014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038. """ if i_density_limit < 1 or i_density_limit > 7: From e8f0ea94d0cea624aa2bd2b32d04fe5e83c30a96 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Thu, 26 Feb 2026 15:04:29 +0000 Subject: [PATCH 17/19] Refactor density limit test to use Physics class method --- tests/unit/test_physics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 9a5d91097..72be4b201 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -2475,7 +2475,7 @@ def test_calculate_density_limit(calculatedensitylimitparam, physics): """ nd_plasma_electron_max_array, nd_plasma_electrons_max = ( - PlasmaDensityLimit().calculate_density_limit( + physics.density_limit.calculate_density_limit( i_density_limit=calculatedensitylimitparam.i_density_limit, b_plasma_toroidal_on_axis=calculatedensitylimitparam.b_plasma_toroidal_on_axis, p_plasma_separatrix_mw=calculatedensitylimitparam.p_plasma_separatrix_mw, From 5dfe58a569e5def8be4e6779b4ccac7c9b656524 Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 2 Mar 2026 13:40:00 +0000 Subject: [PATCH 18/19] Refactor density limit calculation in PlasmaDensityLimit class to use unified method --- process/models/physics/density_limit.py | 94 +++---------------------- 1 file changed, 8 insertions(+), 86 deletions(-) diff --git a/process/models/physics/density_limit.py b/process/models/physics/density_limit.py index a3b0541aa..6bd8611e8 100644 --- a/process/models/physics/density_limit.py +++ b/process/models/physics/density_limit.py @@ -36,99 +36,21 @@ def __init__(self): self.mfile = constants.MFILE def run(self): - nd_plasma_electron_max_array = np.empty((8,)) - - p_perp = ( - physics_variables.p_plasma_separatrix_mw / physics_variables.a_plasma_surface - ) - - # Old ASDEX density limit formula - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - nd_plasma_electron_max_array[0] = self.calculate_asdex_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - q95=physics_variables.q95, - rmajor=physics_variables.rmajor, - prn1=divertor_variables.prn1, - ) - - # Borrass density limit model for ITER (I) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # Borrass et al, ITER-TN-PH-9-6 (1989) - - nd_plasma_electron_max_array[1] = self.calculate_borrass_iter_i_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - q95=physics_variables.q95, - rmajor=physics_variables.rmajor, - prn1=divertor_variables.prn1, - ) - - # Borrass density limit model for ITER (II) - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # This formula is (almost) identical to that in the original routine - # denlim (now deleted). - - nd_plasma_electron_max_array[2] = self.calculate_borrass_iter_ii_density_limit( - p_perp=p_perp, - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - q95=physics_variables.q95, - rmajor=physics_variables.rmajor, - prn1=divertor_variables.prn1, - ) - - # JET edge radiation density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - # qcyl=qstar here, but literature is not clear. - - nd_plasma_electron_max_array[3] = ( - self.calculate_jet_edge_radiation_density_limit( - zeff=physics_variables.n_charge_plasma_effective_vol_avg, - p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, - prn1=divertor_variables.prn1, - qcyl=physics_variables.qstar, - ) - ) - - # JET simplified density limit model - # This applies to the density at the plasma edge, so must be scaled - # to give the density limit applying to the average plasma density. - - nd_plasma_electron_max_array[4] = self.calculate_jet_simple_density_limit( + physics_variables.nd_plasma_electron_max_array, _ = self.calculate_density_limit( b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, + i_density_limit=physics_variables.i_density_limit, p_plasma_separatrix_mw=physics_variables.p_plasma_separatrix_mw, - rmajor=physics_variables.rmajor, + p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, + plasma_current=physics_variables.plasma_current, prn1=divertor_variables.prn1, - ) - - # Hugill-Murakami M.q limit - # qcyl=qstar here, which is okay according to the literature - - nd_plasma_electron_max_array[5] = self.calculate_hugill_murakami_density_limit( - b_plasma_toroidal_on_axis=physics_variables.b_plasma_toroidal_on_axis, - rmajor=physics_variables.rmajor, qcyl=physics_variables.qstar, - ) - - # Greenwald limit - - nd_plasma_electron_max_array[6] = self.calculate_greenwald_density_limit( - c_plasma=physics_variables.plasma_current, rminor=physics_variables.rminor - ) - - nd_plasma_electron_max_array[7] = self.calculate_asdex_new_density_limit( - p_hcd_injected_total_mw=current_drive_variables.p_hcd_injected_total_mw, - c_plasma=physics_variables.plasma_current, q95=physics_variables.q95, - prn1=divertor_variables.prn1, + rmajor=physics_variables.rmajor, + rminor=physics_variables.rminor, + a_plasma_surface=physics_variables.a_plasma_surface, + zeff=physics_variables.n_charge_plasma_effective_vol_avg, ) - physics_variables.nd_plasma_electron_max_array = nd_plasma_electron_max_array - # Calculate beta_norm_max based on i_beta_norm_max try: model = DensityLimitModel(int(physics_variables.i_density_limit)) From 27dd24db3862af443ba349cfa3e1e9f3d4fe932b Mon Sep 17 00:00:00 2001 From: mn3981 Date: Mon, 2 Mar 2026 13:47:55 +0000 Subject: [PATCH 19/19] Refactor import statements in density_limit.py for consistency --- process/models/physics/density_limit.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/process/models/physics/density_limit.py b/process/models/physics/density_limit.py index 6bd8611e8..0819ff76f 100644 --- a/process/models/physics/density_limit.py +++ b/process/models/physics/density_limit.py @@ -3,14 +3,14 @@ import numpy as np -from process import constants -from process import process_output as po +from process.core import constants +from process.core import process_output as po +from process.core.exceptions import ProcessValueError from process.data_structure import ( current_drive_variables, divertor_variables, physics_variables, ) -from process.exceptions import ProcessValueError logger = logging.getLogger(__name__)