feat: Add lambda_q Eich regression 9 (#114)

* Add Eich regression 9

* Backwards compatibility: add reg 9 to named options

* Fix doctest

* Call unitless func for reg 9

* Fix bug: hadn't passed lambda_q to reg 9 in wrapper

* Add tests for lambda_q scalings

* Py3.9 compatibility

cfspopcon/formulas/scrape_off_layer/__init__.py

@@ -1,7 +1,13 @@
 """Routines to calculate the scrape-off-layer conditions and check divertor survivability."""
 
 from .heat_flux_density import calc_B_pol_omp, calc_B_tor_omp, calc_fieldline_pitch_at_omp, calc_parallel_heat_flux_density, calc_q_perp
-from .lambda_q import calc_lambda_q
+from .lambda_q import (
+    calc_lambda_q,
+    calc_lambda_q_with_brunner,
+    calc_lambda_q_with_eich_regression_9,
+    calc_lambda_q_with_eich_regression_14,
+    calc_lambda_q_with_eich_regression_15,
+)
 from .reattachment_models import (
     calc_ionization_volume_from_AUG,
     calc_neutral_flux_density_factor,
@@ -25,6 +31,10 @@
     "two_point_model_fixed_qpart",
     "two_point_model_fixed_tet",
     "calc_lambda_q",
+    "calc_lambda_q_with_brunner",
+    "calc_lambda_q_with_eich_regression_14",
+    "calc_lambda_q_with_eich_regression_15",
+    "calc_lambda_q_with_eich_regression_9",
     "calc_separatrix_electron_density",
     "calc_separatrix_electron_temp",
     "calc_B_pol_omp",

cfspopcon/formulas/scrape_off_layer/lambda_q.py

@@ -15,6 +15,8 @@
         major_radius=ureg.meter,
         B_pol_out_mid=ureg.tesla,
         inverse_aspect_ratio=ureg.dimensionless,
+        magnetic_field_on_axis=ureg.T,
+        q_star=ureg.dimensionless,
         lambda_q_factor=ureg.dimensionless,
     ),
 )
@@ -25,6 +27,8 @@ def calc_lambda_q(
     major_radius: float,
     B_pol_out_mid: float,
     inverse_aspect_ratio: float,
+    magnetic_field_on_axis: float,
+    q_star: float,
     lambda_q_factor: float = 1.0,
 ) -> float:
     """Calculate SOL heat flux decay length (lambda_q) from a scaling.
@@ -36,64 +40,101 @@ def calc_lambda_q(
         major_radius: [m] :term:`glossary link<major_radius>`
         B_pol_out_mid: [T] :term:`glossary link<B_pol_out_mid>`
         inverse_aspect_ratio: [~] :term:`glossary link<inverse_aspect_ratio>`
+        magnetic_field_on_axis: [T] :term:`glossary link<magnetic_field_on_axis>`
+        q_star: [~] :term:`glossary link<q_star>`
         lambda_q_factor: [~] :term:`glossary link<lambda_q_factor>`
 
     Returns:
         :term:`lambda_q` [mm]
     """
     if lambda_q_scaling == LambdaQScaling.Brunner:
-        return lambda_q_factor * float(calc_lambda_q_with_brunner.unitless_func(average_total_pressure))
+        lambda_q = calc_lambda_q_with_brunner.unitless_func(average_total_pressure=average_total_pressure, lambda_q_factor=lambda_q_factor)
+    elif lambda_q_scaling == LambdaQScaling.EichRegression9:
+        lambda_q = calc_lambda_q_with_eich_regression_9.unitless_func(
+            magnetic_field_on_axis=magnetic_field_on_axis,
+            q_star=q_star,
+            power_crossing_separatrix=power_crossing_separatrix,
+            lambda_q_factor=lambda_q_factor,
+        )
     elif lambda_q_scaling == LambdaQScaling.EichRegression14:
-        return lambda_q_factor * float(calc_lambda_q_with_eich_regression_14.unitless_func(B_pol_out_mid))
+        lambda_q = calc_lambda_q_with_eich_regression_14.unitless_func(B_pol_out_mid=B_pol_out_mid, lambda_q_factor=lambda_q_factor)
     elif lambda_q_scaling == LambdaQScaling.EichRegression15:
-        return lambda_q_factor * float(
-            calc_lambda_q_with_eich_regression_15.unitless_func(
-                power_crossing_separatrix, major_radius, B_pol_out_mid, inverse_aspect_ratio
-            )
+        lambda_q = calc_lambda_q_with_eich_regression_15.unitless_func(
+            power_crossing_separatrix=power_crossing_separatrix,
+            major_radius=major_radius,
+            B_pol_out_mid=B_pol_out_mid,
+            inverse_aspect_ratio=inverse_aspect_ratio,
+            lambda_q_factor=lambda_q_factor,
         )
     else:
         raise NotImplementedError(f"No implementation for lambda_q scaling {lambda_q_scaling}")
 
+    return float(lambda_q)
+
 
+@Algorithm.register_algorithm(return_keys=["lambda_q"])
 @wraps_ufunc(
     return_units=dict(lambda_q=ureg.millimeter),
-    input_units=dict(average_total_pressure=ureg.atm),
+    input_units=dict(average_total_pressure=ureg.atm, lambda_q_factor=ureg.dimensionless),
 )
-def calc_lambda_q_with_brunner(average_total_pressure: float) -> float:
+def calc_lambda_q_with_brunner(average_total_pressure: float, lambda_q_factor: float = 1.0) -> float:
     """Return lambda_q according to the Brunner scaling.
 
     Equation 4 in :cite:`brunner_2018_heat_flux`
     """
-    return float(0.91 * average_total_pressure**-0.48)
+    return float(lambda_q_factor * 0.91 * average_total_pressure**-0.48)
 
 
-@wraps_ufunc(return_units=dict(lambda_q=ureg.millimeter), input_units=dict(B_pol_out_mid=ureg.tesla))
-def calc_lambda_q_with_eich_regression_14(B_pol_out_mid: float) -> float:
+@Algorithm.register_algorithm(return_keys=["lambda_q"])
+@wraps_ufunc(
+    return_units=dict(lambda_q=ureg.millimeter),
+    input_units=dict(
+        magnetic_field_on_axis=ureg.T,
+        q_star=ureg.dimensionless,
+        power_crossing_separatrix=ureg.megawatt,
+        lambda_q_factor=ureg.dimensionless,
+    ),
+)
+def calc_lambda_q_with_eich_regression_9(
+    magnetic_field_on_axis: float, q_star: float, power_crossing_separatrix: float, lambda_q_factor: float = 1.0
+) -> float:
+    """Return lambda_q according to Eich regression 9.
+
+    #9 in Table 2 in :cite:`eich_scaling_2013`
+    """
+    return float(lambda_q_factor * 0.7 * magnetic_field_on_axis**-0.77 * q_star**1.05 * power_crossing_separatrix**0.09)
+
+
+@Algorithm.register_algorithm(return_keys=["lambda_q"])
+@wraps_ufunc(return_units=dict(lambda_q=ureg.millimeter), input_units=dict(B_pol_out_mid=ureg.tesla, lambda_q_factor=ureg.dimensionless))
+def calc_lambda_q_with_eich_regression_14(B_pol_out_mid: float, lambda_q_factor: float = 1.0) -> float:
     """Return lambda_q according to Eich regression 14.
 
     #14 in Table 3 in :cite:`eich_scaling_2013`
     """
-    return float(0.63 * B_pol_out_mid**-1.19)
+    return float(lambda_q_factor * 0.63 * B_pol_out_mid**-1.19)
 
 
+@Algorithm.register_algorithm(return_keys=["lambda_q"])
 @wraps_ufunc(
     return_units=dict(lambda_q=ureg.millimeter),
     input_units=dict(
         power_crossing_separatrix=ureg.megawatt,
         major_radius=ureg.meter,
         B_pol_out_mid=ureg.tesla,
         inverse_aspect_ratio=ureg.dimensionless,
+        lambda_q_factor=ureg.dimensionless,
     ),
 )
 def calc_lambda_q_with_eich_regression_15(
-    power_crossing_separatrix: float, major_radius: float, B_pol_out_mid: float, inverse_aspect_ratio: float
+    power_crossing_separatrix: float, major_radius: float, B_pol_out_mid: float, inverse_aspect_ratio: float, lambda_q_factor: float = 1.0
 ) -> float:
     """Return lambda_q according to Eich regression 15.
 
     #15 in Table 3 in :cite:`eich_scaling_2013`
     """
     lambda_q = 1.35 * major_radius**0.04 * B_pol_out_mid**-0.92 * inverse_aspect_ratio**0.42
     if power_crossing_separatrix > 0:
-        return float(lambda_q * power_crossing_separatrix**-0.02)
+        return float(lambda_q_factor * lambda_q * power_crossing_separatrix**-0.02)
     else:
-        return float(lambda_q)
+        return float(lambda_q_factor * lambda_q)

cfspopcon/named_options.py

@@ -63,6 +63,7 @@ class LambdaQScaling(Enum):
     """Options for heat flux decay length scaling."""
 
     Brunner = auto()
+    EichRegression9 = auto()
     EichRegression14 = auto()
     EichRegression15 = auto()
 

cfspopcon/variables.yaml

@@ -153,6 +153,7 @@ average_total_pressure:
   - calc_average_total_pressure
   used_by:
   - calc_lambda_q
+  - calc_lambda_q_with_brunner
 B_pol_out_mid:
   default_units: tesla
   description:
@@ -162,6 +163,8 @@ B_pol_out_mid:
   used_by:
   - calc_fieldline_pitch_at_omp
   - calc_lambda_q
+  - calc_lambda_q_with_eich_regression_14
+  - calc_lambda_q_with_eich_regression_15
   - calc_power_crossing_separatrix_in_electron_channel
 B_t_out_mid:
   default_units: tesla
@@ -740,6 +743,7 @@ inverse_aspect_ratio:
   - calc_plasma_current_from_qstar
   - calc_q_star_from_plasma_current
   - calc_lambda_q
+  - calc_lambda_q_with_eich_regression_15
 invmu_0_dLedR:
   default_units: dimensionless
   description:
@@ -838,6 +842,10 @@ lambda_q:
   - The near-scrape-off-layer heat-flux-density decay length.
   set_by:
   - calc_lambda_q
+  - calc_lambda_q_with_brunner
+  - calc_lambda_q_with_eich_regression_9
+  - calc_lambda_q_with_eich_regression_14
+  - calc_lambda_q_with_eich_regression_15
   used_by:
   - calc_parallel_heat_flux_density
   - calc_q_perp
@@ -850,6 +858,10 @@ lambda_q_factor:
   set_by: []
   used_by:
   - calc_lambda_q
+  - calc_lambda_q_with_brunner
+  - calc_lambda_q_with_eich_regression_9
+  - calc_lambda_q_with_eich_regression_14
+  - calc_lambda_q_with_eich_regression_15
 lambda_q_scaling:
   default_units: null
   description:
@@ -897,6 +909,8 @@ magnetic_field_on_axis:
   - calc_beta_normalized
   - calc_troyon_limit
   - calc_B_tor_omp
+  - calc_lambda_q
+  - calc_lambda_q_with_eich_regression_9
   - calc_SepOS_L_mode_density_limit
   - calc_SepOS_LH_transition
   - calc_SepOS_ideal_MHD_limit
@@ -940,6 +954,7 @@ major_radius:
   - calc_parallel_heat_flux_density
   - calc_q_perp
   - calc_lambda_q
+  - calc_lambda_q_with_eich_regression_15
   - calc_SepOS_L_mode_density_limit
   - calc_SepOS_LH_transition
   - calc_SepOS_ideal_MHD_limit
@@ -1406,6 +1421,8 @@ power_crossing_separatrix:
   - calc_parallel_heat_flux_density
   - calc_q_perp
   - calc_lambda_q
+  - calc_lambda_q_with_eich_regression_9
+  - calc_lambda_q_with_eich_regression_15
   - calc_ratio_P_LH
   - calc_ratio_P_LI
 pumping_duct_neutral_pressure:
@@ -1456,6 +1473,8 @@ q_star:
   - calc_PBpRnSq
   - calc_bootstrap_fraction
   - calc_plasma_current_from_qstar
+  - calc_lambda_q
+  - calc_lambda_q_with_eich_regression_9
 radas_dir:
   default_units: null
   description:

docs/conf.py

@@ -51,6 +51,9 @@
     r"https://www.ipp.mpg.de/16701/jet",
     r"https://iopscience.iop.org/article/10.1088/1009-0630/13/1/01",
     r"https://www-internal.psfc.mit.edu/research/alcator/data/fst_cmod.pdf",
+    # These bib resources fail due to "403 Client Error: Forbidden for url"
+    r"https://doi.org/10.1103/PhysRevLett.121.055001",
+    r"https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.121.055001",
 ]
 linkcheck_retries = 5
 linkcheck_timeout = 120

tests/test_scrape_off_layer_model.py

@@ -0,0 +1,124 @@
+import numpy as np
+import pytest
+
+from cfspopcon import formulas
+from cfspopcon.named_options import LambdaQScaling
+from cfspopcon.unit_handling import ureg
+
+lambda_q_tests = {
+    LambdaQScaling.Brunner: 0.4332283874128845 * ureg.mm,
+    LambdaQScaling.EichRegression14: 0.20533809707365488 * ureg.mm,
+    LambdaQScaling.EichRegression15: 0.34842497310813536 * ureg.mm,
+    LambdaQScaling.EichRegression9: 0.5865460692254366 * ureg.mm,
+}
+
+
+@pytest.fixture()
+def average_total_pressure():
+    return 732028.9793 * ureg.Pa
+
+
+@pytest.fixture()
+def power_crossing_separatrix():
+    return 25.57417052 * ureg.MW
+
+
+@pytest.fixture()
+def major_radius():
+    return 1.85 * ureg.m
+
+
+@pytest.fixture()
+def B_pol_out_mid():
+    return 3.052711915 * ureg.T
+
+
+@pytest.fixture()
+def inverse_aspect_ratio():
+    return 0.3081000000
+
+
+@pytest.fixture()
+def magnetic_field_on_axis():
+    return 12.20000000 * ureg.T
+
+
+@pytest.fixture()
+def q_star():
+    return 3.290275716
+
+
+@pytest.fixture()
+def lambda_q_factor():
+    return 1.23
+
+
+@pytest.mark.parametrize(["scaling", "result"], lambda_q_tests.items(), ids=[key.name for key in lambda_q_tests.keys()])
+def test_lambda_q_scalings(
+    scaling,
+    result,
+    average_total_pressure,
+    power_crossing_separatrix,
+    major_radius,
+    B_pol_out_mid,
+    inverse_aspect_ratio,
+    magnetic_field_on_axis,
+    q_star,
+    lambda_q_factor,
+):
+    lambda_q = formulas.scrape_off_layer.calc_lambda_q(
+        lambda_q_scaling=scaling,
+        average_total_pressure=average_total_pressure,
+        power_crossing_separatrix=power_crossing_separatrix,
+        major_radius=major_radius,
+        B_pol_out_mid=B_pol_out_mid,
+        inverse_aspect_ratio=inverse_aspect_ratio,
+        magnetic_field_on_axis=magnetic_field_on_axis,
+        q_star=q_star,
+        lambda_q_factor=lambda_q_factor,
+    )
+
+    assert np.isclose(lambda_q, result)
+
+
+@pytest.mark.parametrize(["scaling", "result"], lambda_q_tests.items(), ids=[key.name for key in lambda_q_tests.keys()])
+def test_lambda_q_scalings_with_algorithms(
+    scaling,
+    result,
+    average_total_pressure,
+    power_crossing_separatrix,
+    major_radius,
+    B_pol_out_mid,
+    inverse_aspect_ratio,
+    magnetic_field_on_axis,
+    q_star,
+    lambda_q_factor,
+):
+    if scaling == LambdaQScaling.Brunner:
+        lambda_q = formulas.scrape_off_layer.lambda_q.calc_lambda_q_with_brunner(
+            average_total_pressure=average_total_pressure, lambda_q_factor=lambda_q_factor
+        )
+    elif scaling == LambdaQScaling.EichRegression14:
+        lambda_q = formulas.scrape_off_layer.lambda_q.calc_lambda_q_with_eich_regression_14(
+            B_pol_out_mid=B_pol_out_mid,
+            lambda_q_factor=lambda_q_factor,
+        )
+    elif scaling == LambdaQScaling.EichRegression15:
+        lambda_q = formulas.scrape_off_layer.lambda_q.calc_lambda_q_with_eich_regression_15(
+            power_crossing_separatrix=power_crossing_separatrix,
+            major_radius=major_radius,
+            B_pol_out_mid=B_pol_out_mid,
+            inverse_aspect_ratio=inverse_aspect_ratio,
+            lambda_q_factor=lambda_q_factor,
+        )
+    elif scaling == LambdaQScaling.EichRegression9:
+        lambda_q = formulas.scrape_off_layer.lambda_q.calc_lambda_q_with_eich_regression_9(
+            magnetic_field_on_axis=magnetic_field_on_axis,
+            q_star=q_star,
+            power_crossing_separatrix=power_crossing_separatrix,
+            lambda_q_factor=lambda_q_factor,
+        )
+    else:
+        raise NotImplementedError(f"Add the algorithm for {scaling.name}.")
+
+    assert np.isclose(lambda_q, result)