I mode power threshold updates (#58)

* Transfer PLI scaling to new branch * Formatting * Add datatype * Remove radas * Handle confinement_power_scaling as a named option * added confinement_power_scaling to physics glossary --------- Co-authored-by: Tom Body <[email protected]> Co-authored-by: aarooshgr <[email protected]>
cfs-energy · Jun 24, 2024 · 887bd82 · 887bd82
1 parent 9040a0c
commit 887bd82
Show file tree

Hide file tree

Showing 6 changed files with 135 additions and 17 deletions.
diff --git a/cfspopcon/default_units.yaml b/cfspopcon/default_units.yaml
@@ -12,6 +12,7 @@ beta_poloidal: ""
 beta_toroidal: ""
 beta: ""
 bootstrap_fraction: ""
+confinement_power_scaling: null
 confinement_threshold_scalar: ""
 confinement_time_scalar: ""
 core_radiated_power_fraction: ""
@@ -66,8 +67,8 @@ lambda_q: "mm"
 lengyel_overestimation_factor: ""
 loop_voltage: "V"
 magnetic_field_on_axis: "T"
-mean_ion_charge_state: ""
 major_radius: "m"
+mean_ion_charge_state: ""
 minimum_core_radiated_fraction: ""
 minor_radius: "m"
 neoclassical_loop_resistivity: "m * ohm"
@@ -84,6 +85,7 @@ P_fusion: "MW"
 P_in: "MW"
 P_launched: "MW"
 P_LH_thresh: "MW"
+P_LI_thresh: "MW"
 P_neutron: "MW"
 P_ohmic: "MW"
 P_radiated_by_core_radiator: "MW"
@@ -113,6 +115,7 @@ radas_dir: null
 radiated_power_method: null
 radiated_power_scalar: ""
 ratio_of_P_SOL_to_P_LH: ""
+ratio_of_P_SOL_to_P_LI: ""
 rho_s_pol: "mm"
 rho_star: ""
 rho: ""
@@ -129,6 +132,8 @@ SOL_power_loss_fraction: ""
 spitzer_resistivity: "m * ohm"
 summed_impurity_density: "n19"
 surface_area: "m**2"
+sustainment_power_in_electron_channel: "MW"
+sustainment_power_in_ion_channel: "MW"
 target_electron_density: "n19"
 target_electron_flux: "m**-2 s**-1"
 target_electron_temp: "eV"
@@ -145,6 +150,4 @@ upstream_electron_density: "n19"
 upstream_electron_temp: "eV"
 vertical_minor_radius: "m"
 z_effective: ""
-zeff_change_from_core_rad: ""
-sustainment_power_in_electron_channel: "MW"
-sustainment_power_in_ion_channel: "MW"
+zeff_change_from_core_rad: ""
diff --git a/cfspopcon/formulas/separatrix_conditions/threshold_power.py b/cfspopcon/formulas/separatrix_conditions/threshold_power.py
@@ -1,6 +1,7 @@
 """Calculate the threshold values for the power crossing the separatrix."""
 
 from ...algorithm_class import Algorithm
+from ...named_options import ConfinementPowerScaling
 from ...unit_handling import ureg, wraps_ufunc
 
 
@@ -15,6 +16,7 @@
         surface_area=ureg.m**2,
         fuel_average_mass_number=ureg.amu,
         average_electron_density=ureg.n19,
+        confinement_power_scaling=None,
         confinement_threshold_scalar=ureg.dimensionless,
     ),
 )
@@ -26,6 +28,7 @@ def calc_LH_transition_threshold_power(
     surface_area: float,
     fuel_average_mass_number: float,
     average_electron_density: float,
+    confinement_power_scaling: ConfinementPowerScaling = ConfinementPowerScaling.H_mode_Martin,
     confinement_threshold_scalar: float = 1.0,
 ) -> float:
     """Calculate the threshold power (crossing the separatrix) to transition into H-mode.
@@ -41,6 +44,7 @@ def calc_LH_transition_threshold_power(
         surface_area: [m^2] :term:`glossary link<surface_area>`
         fuel_average_mass_number: [amu] :term:`glossary link<fuel_average_mass_number>`
         average_electron_density: [1e19 m^-3] :term:`glossary link<average_electron_density>`
+        confinement_power_scaling: [~] :term:`glossary link<confinement_power_scaling>`
         confinement_threshold_scalar: [~] :term:`glossary link<confinement_threshold_scalar>`
 
     Returns:
@@ -54,16 +58,27 @@ def _calc_Martin_LH_threshold(electron_density: float) -> float:
             _DEUTERIUM_MASS_NUMBER / fuel_average_mass_number
         )
 
-    neMin19 = (
-        0.7 * (plasma_current**0.34) * (magnetic_field_on_axis**0.62) * (minor_radius**-0.95) * ((major_radius / minor_radius) ** 0.4)
-    )
+    if confinement_power_scaling == ConfinementPowerScaling.H_mode_Martin:
+
+        neMin19 = (
+            0.7
+            * (plasma_current**0.34)
+            * (magnetic_field_on_axis**0.62)
+            * (minor_radius**-0.95)
+            * ((major_radius / minor_radius) ** 0.4)
+        )
+
+        if average_electron_density < neMin19:
+            P_LH_thresh = _calc_Martin_LH_threshold(electron_density=neMin19)
+            return float(P_LH_thresh * (neMin19 / average_electron_density) ** 2.0) * confinement_threshold_scalar
+        else:
+            P_LH_thresh = _calc_Martin_LH_threshold(electron_density=average_electron_density)
+            return P_LH_thresh * confinement_threshold_scalar
 
-    if average_electron_density < neMin19:
-        P_LH_thresh = _calc_Martin_LH_threshold(electron_density=neMin19)
-        return float(P_LH_thresh * (neMin19 / average_electron_density) ** 2.0) * confinement_threshold_scalar
     else:
-        P_LH_thresh = _calc_Martin_LH_threshold(electron_density=average_electron_density)
-        return P_LH_thresh * confinement_threshold_scalar
+        raise NotImplementedError(
+            f"No implementation for calc_LH_transition_threshold_power with confinement_power_scaling = {confinement_power_scaling} (type ({type(confinement_power_scaling)}))"
+        )
 
 
 calc_ratio_P_LH = Algorithm.from_single_function(
@@ -74,24 +89,79 @@ def _calc_Martin_LH_threshold(electron_density: float) -> float:
 @Algorithm.register_algorithm(return_keys=["P_LI_thresh"])
 @wraps_ufunc(
     return_units=dict(P_LI_thresh=ureg.MW),
-    input_units=dict(plasma_current=ureg.MA, average_electron_density=ureg.n19, confinement_threshold_scalar=ureg.dimensionless),
+    input_units=dict(
+        plasma_current=ureg.MA,
+        magnetic_field_on_axis=ureg.T,
+        surface_area=ureg.m**2,
+        average_electron_density=ureg.n19,
+        confinement_power_scaling=None,
+        confinement_threshold_scalar=ureg.dimensionless,
+    ),
 )
 def calc_LI_transition_threshold_power(
-    plasma_current: float, average_electron_density: float, confinement_threshold_scalar: float = 1.0
+    plasma_current: float,
+    magnetic_field_on_axis: float,
+    surface_area: float,
+    average_electron_density: float,
+    confinement_power_scaling: ConfinementPowerScaling = ConfinementPowerScaling.I_mode_HubbardNF17,
+    confinement_threshold_scalar: float = 1.0,
 ) -> float:
     """Calculate the threshold power (crossing the separatrix) to transition into I-mode.
 
-    Note: uses scaling described in Fig 5 of ref :cite:`hubbard_threshold_2012`
+    Note:
+        "AUG" option is inspired from :cite:`ryter_i-mode_2016` and :cite:`Happel_2017`
+        "HubbardNF17" uses scaling described in Fig 6 of :cite:`hubbard_threshold_2017`
+        "HubbardNF12" uses scaling described in Fig 5 of ref :cite:`hubbard_threshold_2012`
 
     Args:
         plasma_current: [MA] :term:`glossary link<plasma_current>`
+        magnetic_field_on_axis: [T] :term:`glossary link<magnetic_field_on_axis>`
+        surface_area: [m^2] :term:`glossary link<surface_area>`
         average_electron_density: [1e19 m^-3] :term:`glossary link<average_electron_density>`
+        confinement_power_scaling: [~] :term:`glossary link<confinement_power_scaling>`
         confinement_threshold_scalar: [~] :term:`glossary link<confinement_threshold_scalar>`
 
     Returns:
         :term:`P_LI_thresh` [MW]
     """
-    return float(2.11 * plasma_current**0.94 * ((average_electron_density / 10.0) ** 0.65)) * confinement_threshold_scalar
+
+    def _calc_AUG_LI_threshold(
+        average_electron_density: float, magnetic_field_on_axis: float, surface_area: float, confinement_threshold_scalar: float
+    ) -> float:
+        return (
+            float(0.14 * (average_electron_density / 10) * (magnetic_field_on_axis / 2.4) ** 0.39 * surface_area)
+            * confinement_threshold_scalar
+        )
+
+    def _calc_HubbardNF17_LI_threshold(
+        average_electron_density: float, magnetic_field_on_axis: float, surface_area: float, confinement_threshold_scalar: float
+    ) -> float:
+        return (
+            float(0.162 * (average_electron_density / 10) * (magnetic_field_on_axis**0.262) * surface_area) * confinement_threshold_scalar
+        )
+
+    def _calc_HubbardNF12_LI_threshold(
+        average_electron_density: float, plasma_current: float, confinement_threshold_scalar: float
+    ) -> float:
+        return float(2.11 * plasma_current**0.94 * ((average_electron_density / 10.0) ** 0.65)) * confinement_threshold_scalar
+
+    if confinement_power_scaling == ConfinementPowerScaling.I_mode_AUG:
+        P_LI_thresh = _calc_AUG_LI_threshold(average_electron_density, magnetic_field_on_axis, surface_area, confinement_threshold_scalar)
+
+    elif confinement_power_scaling == ConfinementPowerScaling.I_mode_HubbardNF17:
+        P_LI_thresh = _calc_HubbardNF17_LI_threshold(
+            average_electron_density, magnetic_field_on_axis, surface_area, confinement_threshold_scalar
+        )
+
+    elif confinement_power_scaling == ConfinementPowerScaling.I_mode_HubbardNF12:
+        P_LI_thresh = _calc_HubbardNF12_LI_threshold(average_electron_density, plasma_current, confinement_threshold_scalar)
+
+    else:
+        raise NotImplementedError(
+            f"No implementation for calc_LI_transition_threshold_power with confinement_power_scaling = {confinement_power_scaling} (type ({type(confinement_power_scaling)}))"
+        )
+
+    return float(P_LI_thresh)
 
 
 calc_ratio_P_LI = Algorithm.from_single_function(

diff --git a/cfspopcon/helpers.py b/cfspopcon/helpers.py
@@ -5,6 +5,7 @@
 
 from .named_options import (
     AtomicSpecies,
+    ConfinementPowerScaling,
     LambdaQScaling,
     MomentumLossFunction,
     ProfileForm,
@@ -30,6 +31,8 @@ def convert_named_options(key: str, val: Any) -> Any:  # noqa: PLR0911
         return MomentumLossFunction[val]
     elif key == "radiation_method":
         return RadiationMethod[val]
+    elif key == "confinement_power_scaling":
+        return ConfinementPowerScaling[val]
     else:
         # If the key doesn't match, don't convert the value
         return val

diff --git a/cfspopcon/named_options.py b/cfspopcon/named_options.py
@@ -64,3 +64,12 @@ class LambdaQScaling(Enum):
     Brunner = auto()
     EichRegression14 = auto()
     EichRegression15 = auto()
+
+
+class ConfinementPowerScaling(Enum):
+    """Options for which confinement threshold power scaling to use."""
+
+    I_mode_AUG = auto()
+    I_mode_HubbardNF17 = auto()
+    I_mode_HubbardNF12 = auto()
+    H_mode_Martin = auto()
diff --git a/docs/doc_sources/physics_glossary.rst b/docs/doc_sources/physics_glossary.rst
@@ -454,4 +454,7 @@ Physics Glossary
     The power in the ion channel required to maintain the ion temperature gradient at the separatrix.
 
   sustainment_power_in_electron_channel
-    The power in the electron channel required to maintain the electron temperature at the separatrix.
+    The power in the electron channel required to maintain the electron temperature at the separatrix.
+
+  confinement_power_scaling
+    A :class:`~cfspopcon.named_options.ConfinementPowerScaling` indicating which confinement mode power-threshold scaling to use.
diff --git a/docs/refs.bib b/docs/refs.bib
@@ -288,6 +288,36 @@ @article{hubbard_threshold_2012
 	pages = {114009},
 }
 
+@article{hubbard_threshold_2017,
+	doi = {10.1088/1741-4326/aa8570},
+	url = {https://dx.doi.org/10.1088/1741-4326/aa8570},
+	year = {2017},
+	month = {oct},
+	publisher = {IOP Publishing},
+	volume = {57},
+	number = {12},
+	pages = {126039},
+	author = {A.E. Hubbard and S.-G. Baek and D. Brunner and A.J. Creely and I. Cziegler and E. Edlund and J.W. Hughes and B. LaBombard and Y. Lin and Z. Liu and E.S. Marmar and M.L. Reinke and J.E. Rice and B. Sorbom and C. Sung and J. Terry and C. Theiler and E.A. Tolman and J.R. Walk and A.E. White and D. Whyte and S.M. Wolfe and S. Wukitch and X.Q. Xu and the Alcator C-Mod team},
+	title = {Physics and performance of the I-mode regime over an expanded operating space on Alcator C-Mod},
+	journal = {Nuclear Fusion},
+	abstract = {New results on the I-mode regime of operation on the Alcator C-Mod tokamak are reported. This ELM-free regime features high energy confinement and a steep temperature pedestal, while particle confinement remains at L-mode levels, giving stationary density and avoiding impurity accumulation. I-mode has now been obtained over nearly all of the magnetic fields and currents possible in this high field tokamak (Ip 0.55–1.7 MA, BT 2.8–8 T) using a configuration with B  ×  ∇B drift away from the X-point. Results at 8 T confirm that the L–I power threshold varies only weakly with BT, and that the power range for I-mode increases with BT; no 8 T discharges transitioned to H-mode. Parameter dependences of energy confinement are investigated. Core transport simulations are giving insight into the observed turbulence reduction, profile stiffness and confinement improvement. Pedestal models explain the observed stability to ELMs, and can simulate the observed weakly coherent mode. Conditions for I–H transitions have complex dependences on density as well as power. I-modes have now been maintained in near-DN configurations, leading to improved divertor power flux sharing. Prospects for I-mode on future fusion devices such as ITER and ARC are encouraging. Further experiments on other tokamaks are needed to improve confidence in extrapolation.}
+}
+
+@article{Happel_2017,
+	doi = {10.1088/0741-3335/59/1/014004},
+	url = {https://dx.doi.org/10.1088/0741-3335/59/1/014004},
+	year = {2016},
+	month = {oct},
+	publisher = {IOP Publishing},
+	volume = {59},
+	number = {1},
+	pages = {014004},
+	author = {T Happel and P Manz and F Ryter and M Bernert and M Dunne and P Hennequin and A Hetzenecker and U Stroth and G D Conway and L Guimarais and C Honoré and E Viezzer and The ASDEX Upgrade Team},
+	title = {The I-mode confinement regime at ASDEX Upgrade: global properties and characterization of strongly intermittent density fluctuations},
+	journal = {Plasma Physics and Controlled Fusion},
+	abstract = {Properties of the I-mode confinement regime on the ASDEX Upgrade tokamak are summarized. A weak dependence of the power threshold for the L-I transition on the toroidal magnetic field strength is found. During improved confinement, the edge radial electric field well deepens. Stability calculations show that the I-mode pedestal is peeling-ballooning stable. Turbulence investigations reveal strongly intermittent density fluctuations linked to the weakly coherent mode in the confined plasma, which become stronger as the confinement quality increases. Across all investigated structure sizes (–12 cm−1, with  the perpendicular wavenumber of turbulent density fluctuations), the intermittent turbulence bursts are observed. Comparison with bolometry data shows that they move poloidally toward the X-point and finally end up in the divertor. This might be indicative that they play a role in inhibiting the density profile growth, such that no pedestal is formed in the edge density profile.}
+}
+
 @article{angioni_scaling_2007,
 	title = {Scaling of density peaking in {H}-mode plasmas based on a combined database of {AUG} and {JET} observations},
 	volume = {47},