L/H-mode confinement time consistent with fLH (#70)

* Add switch for L-mode confinement below PLH

* Add test for switch to L-mode below LH power threshold

* ruff ruff

cfspopcon/formulas/energy_confinement/__init__.py

@@ -3,12 +3,16 @@
 from .plasma_stored_energy import calc_plasma_stored_energy
 from .read_energy_confinement_scalings import ConfinementScaling, read_confinement_scalings
 from .solve_for_input_power import solve_energy_confinement_scaling_for_input_power
-from .switch_confinement_scaling_on_threshold import switch_to_linearised_ohmic_confinement_below_threshold
+from .switch_confinement_scaling_on_threshold import (
+    switch_to_L_mode_confinement_below_threshold,
+    switch_to_linearised_ohmic_confinement_below_threshold,
+)
 
 __all__ = [
     "calc_plasma_stored_energy",
     "read_confinement_scalings",
     "ConfinementScaling",
     "solve_energy_confinement_scaling_for_input_power",
     "switch_to_linearised_ohmic_confinement_below_threshold",
+    "switch_to_L_mode_confinement_below_threshold",
 ]

cfspopcon/formulas/energy_confinement/switch_confinement_scaling_on_threshold.py

@@ -72,4 +72,69 @@ def switch_to_linearised_ohmic_confinement_below_threshold(
     return (energy_confinement_time, P_in, SOC_LOC_ratio)
 
 
-# TODO implement switch to L-mode below PLH
+@Algorithm.register_algorithm(return_keys=["energy_confinement_time", "P_in"])
+def switch_to_L_mode_confinement_below_threshold(
+    plasma_stored_energy: Unitfull,
+    energy_confinement_time: Unitfull,
+    P_in: Unitfull,
+    average_electron_density: Unitfull,
+    confinement_time_scalar: Unitfull,
+    plasma_current: Unitfull,
+    magnetic_field_on_axis: Unitfull,
+    major_radius: Unitfull,
+    areal_elongation: Unitfull,
+    separatrix_elongation: Unitfull,
+    inverse_aspect_ratio: Unitfull,
+    average_ion_mass: Unitfull,
+    triangularity_psi95: Unitfull,
+    separatrix_triangularity: Unitfull,
+    q_star: Unitfull,
+    ratio_of_P_SOL_to_P_LH: Unitfull,
+    energy_confinement_scaling_for_L_mode: Unitfull,
+) -> tuple[Unitfull, ...]:
+    """Switch to the L-mode scaling if Psol < PLH.
+
+    Args:
+        plasma_stored_energy: :term:`glossary link<plasma_stored_energy>`
+        energy_confinement_time: :term:`glossary link<energy_confinement_time>`
+        P_in: :term:`glossary link<P_in>`
+        average_electron_density: :term:`glossary link<average_electron_density>`
+        confinement_time_scalar: :term:`glossary link<confinement_time_scalar>`
+        plasma_current: :term:`glossary link<plasma_current>`
+        magnetic_field_on_axis: :term:`glossary link<magnetic_field_on_axis>`
+        major_radius: :term:`glossary link<major_radius>`
+        areal_elongation: :term:`glossary link<areal_elongation>`
+        separatrix_elongation: :term:`glossary link<separatrix_elongation>`
+        inverse_aspect_ratio: :term:`glossary link<inverse_aspect_ratio>`
+        average_ion_mass: :term:`glossary link<average_ion_mass>`
+        triangularity_psi95: :term:`glossary link<triangularity_psi95>`
+        separatrix_triangularity: :term:`glossary link<separatrix_triangularity>`
+        q_star: :term:`glossary link<q_star>`
+        ratio_of_P_SOL_to_P_LH: :term:`glossary link<ratio_of_P_SOL_to_P_LH>`
+        energy_confinement_scaling_for_L_mode: :term:`glossary link<energy_confinement_scaling_for_L_mode>`
+
+    Returns:
+        :term:`energy_confinement_time`, :term:`P_in`
+    """
+    tau_e_L_mode, P_in_L_mode = solve_energy_confinement_scaling_for_input_power(
+        confinement_time_scalar=confinement_time_scalar,
+        plasma_current=plasma_current,
+        magnetic_field_on_axis=magnetic_field_on_axis,
+        average_electron_density=average_electron_density,
+        major_radius=major_radius,
+        areal_elongation=areal_elongation,
+        separatrix_elongation=separatrix_elongation,
+        inverse_aspect_ratio=inverse_aspect_ratio,
+        average_ion_mass=average_ion_mass,
+        triangularity_psi95=triangularity_psi95,
+        separatrix_triangularity=separatrix_triangularity,
+        plasma_stored_energy=plasma_stored_energy,
+        q_star=q_star,
+        energy_confinement_scaling=energy_confinement_scaling_for_L_mode,
+    )
+
+    # Use L-mode confinement if Psol < PLH
+    energy_confinement_time = xr.where(ratio_of_P_SOL_to_P_LH < 1.0, tau_e_L_mode, energy_confinement_time)  # type:ignore[no-untyped-call]
+    P_in = xr.where(ratio_of_P_SOL_to_P_LH < 1.0, P_in_L_mode, P_in)  # type:ignore[no-untyped-call]
+
+    return energy_confinement_time, P_in

cfspopcon/variables.yaml

@@ -20,6 +20,7 @@ areal_elongation:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_plasma_volume
   - calc_plasma_surface_area
   - calc_plasma_poloidal_circumference
@@ -58,6 +59,7 @@ average_electron_density:
   - calc_plasma_stored_energy
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_impurity_charge_state
   - calc_zeff_and_dilution_due_to_impurities
   - calc_normalised_collisionality
@@ -114,6 +116,7 @@ average_ion_mass:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_alpha_t
   - calc_edge_collisionality
   - calc_rho_star
@@ -257,6 +260,7 @@ confinement_time_scalar:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
 core_impurity_species:
   default_units: null
   description:
@@ -469,6 +473,14 @@ energy_confinement_scaling:
   set_by: []
   used_by:
   - solve_energy_confinement_scaling_for_input_power
+energy_confinement_scaling_for_L_mode:
+  default_units: null
+  description:
+  - Which :math:`\tau_e` energy confinement scaling should be used if :math:`P_{SOL}
+    < P_{LH}`. Should match an L-mode confinement scaling in `cfspopcon.formulas.energy_confinement::energy_confinement_scalings.yaml`.
+  set_by: []
+  used_by:
+  - switch_to_L_mode_confinement_below_threshold
 energy_confinement_time:
   default_units: second
   description:
@@ -477,8 +489,10 @@ energy_confinement_time:
   set_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   used_by:
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_triple_product
 external_flux:
   default_units: weber
@@ -708,6 +722,7 @@ inverse_aspect_ratio:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_plasma_volume
   - calc_plasma_surface_area
   - calc_minor_radius_from_inverse_aspect_ratio
@@ -841,6 +856,7 @@ magnetic_field_on_axis:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_synchrotron_radiation
   - calc_intrinsic_radiated_power_from_core
   - calc_PB_over_R
@@ -866,6 +882,7 @@ major_radius:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_plasma_volume
   - calc_plasma_surface_area
   - calc_minor_radius_from_inverse_aspect_ratio
@@ -1075,8 +1092,10 @@ P_in:
   set_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   used_by:
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_fusion_gain
   - calc_f_rad_core
   - require_P_rad_less_than_P_in
@@ -1260,6 +1279,7 @@ plasma_current:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_greenwald_density_limit
   - calc_inductive_plasma_current
   - calc_internal_flux
@@ -1284,6 +1304,7 @@ plasma_stored_energy:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
 plasma_volume:
   default_units: meter ** 3
   description:
@@ -1380,6 +1401,7 @@ q_star:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_normalised_collisionality
   - calc_PBpRnSq
   - calc_bootstrap_fraction
@@ -1417,7 +1439,8 @@ ratio_of_P_SOL_to_P_LH:
     power.
   set_by:
   - calc_ratio_P_LH
-  used_by: []
+  used_by:
+  - switch_to_L_mode_confinement_below_threshold
 ratio_of_P_SOL_to_P_LI:
   default_units: dimensionless
   description:
@@ -1541,6 +1564,7 @@ separatrix_elongation:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_vertical_minor_radius_from_elongation_and_minor_radius
   - calc_synchrotron_radiation
   - calc_intrinsic_radiated_power_from_core
@@ -1553,6 +1577,7 @@ separatrix_triangularity:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
 SepOS_density_limit:
   default_units: dimensionless
   description:
@@ -1763,6 +1788,7 @@ triangularity_psi95:
   used_by:
   - solve_energy_confinement_scaling_for_input_power
   - switch_to_linearised_ohmic_confinement_below_threshold
+  - switch_to_L_mode_confinement_below_threshold
   - calc_separatrix_triangularity_from_triangularity95
   - calc_f_shaping_for_qstar
   - calc_cylindrical_edge_safety_factor

docs/doc_sources/physics_glossary.rst

@@ -137,6 +137,9 @@ Physics Glossary
   energy_confinement_scaling
     Which :math:`\tau_e` energy confinement scaling should be used. Should match a confinement scaling in `cfspopcon.formulas.energy_confinement::energy_confinement_scalings.yaml`.
 
+  energy_confinement_scaling_for_L_mode
+    Which :math:`\tau_e` energy confinement scaling should be used if :math:`P_{SOL} < P_{LH}`. Should match an L-mode confinement scaling in `cfspopcon.formulas.energy_confinement::energy_confinement_scalings.yaml`.
+
   energy_confinement_time
     A characteristic time which gives the rate at which the plasma loses energy. In steady-state, :math:`\tau_e=W_p / P_in`.
 

tests/test_confinement_switch.py

@@ -0,0 +1,54 @@
+from cfspopcon.formulas.energy_confinement.switch_confinement_scaling_on_threshold import switch_to_L_mode_confinement_below_threshold
+from cfspopcon.formulas.energy_confinement.solve_for_input_power import solve_energy_confinement_scaling_for_input_power
+from cfspopcon.formulas.energy_confinement.read_energy_confinement_scalings import read_confinement_scalings
+from cfspopcon.unit_handling import ureg, magnitude_in_units
+import numpy as np
+
+
+def test_switch_to_L_mode_confinement_below_threshold():
+    kwargs = dict(
+        plasma_stored_energy=20.0 * ureg.MJ,
+        average_electron_density=25.0 * ureg.n19,
+        confinement_time_scalar=1.0,
+        plasma_current=8.7 * ureg.MA,
+        magnetic_field_on_axis=12.16 * ureg.T,
+        major_radius=1.85 * ureg.m,
+        areal_elongation=1.75,
+        separatrix_elongation=1.96,
+        inverse_aspect_ratio=0.308,
+        average_ion_mass=2.5 * ureg.amu,
+        triangularity_psi95=0.3,
+        separatrix_triangularity=0.54,
+        q_star=3.29,
+    )
+
+    read_confinement_scalings()
+
+    tau_E_H_mode, P_in_H_mode = solve_energy_confinement_scaling_for_input_power(
+        **kwargs,
+        energy_confinement_scaling="ITER98y2",
+    )
+
+    tau_E_L_mode, _ = solve_energy_confinement_scaling_for_input_power(
+        **kwargs,
+        energy_confinement_scaling="ITER89P",
+    )
+
+    tau_E_should_be_L_mode, _ = switch_to_L_mode_confinement_below_threshold(
+        **kwargs,
+        ratio_of_P_SOL_to_P_LH=0.90,
+        energy_confinement_time=tau_E_H_mode,
+        P_in=P_in_H_mode,
+        energy_confinement_scaling_for_L_mode="ITER89P",
+    )
+
+    tau_E_should_be_H_mode, _ = switch_to_L_mode_confinement_below_threshold(
+        **kwargs,
+        ratio_of_P_SOL_to_P_LH=1.2,
+        energy_confinement_time=tau_E_H_mode,
+        P_in=P_in_H_mode,
+        energy_confinement_scaling_for_L_mode="ITER89P",
+    )
+
+    np.testing.assert_allclose(magnitude_in_units(tau_E_should_be_L_mode, ureg.s), magnitude_in_units(tau_E_L_mode, ureg.s), rtol=1e-3)
+    np.testing.assert_allclose(magnitude_in_units(tau_E_should_be_H_mode, ureg.s), magnitude_in_units(tau_E_H_mode, ureg.s), rtol=1e-3)