diff --git a/cfspopcon/formulas/scrape_off_layer/__init__.py b/cfspopcon/formulas/scrape_off_layer/__init__.py index 1419d67..b304429 100644 --- a/cfspopcon/formulas/scrape_off_layer/__init__.py +++ b/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", diff --git a/cfspopcon/formulas/scrape_off_layer/lambda_q.py b/cfspopcon/formulas/scrape_off_layer/lambda_q.py index a138bc7..ad749f0 100644 --- a/cfspopcon/formulas/scrape_off_layer/lambda_q.py +++ b/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,46 +40,82 @@ def calc_lambda_q( major_radius: [m] :term:`glossary link` B_pol_out_mid: [T] :term:`glossary link` inverse_aspect_ratio: [~] :term:`glossary link` + magnetic_field_on_axis: [T] :term:`glossary link` + q_star: [~] :term:`glossary link` lambda_q_factor: [~] :term:`glossary link` 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( @@ -83,10 +123,11 @@ def calc_lambda_q_with_eich_regression_14(B_pol_out_mid: float) -> float: 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. @@ -94,6 +135,6 @@ def calc_lambda_q_with_eich_regression_15( """ 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) diff --git a/cfspopcon/named_options.py b/cfspopcon/named_options.py index 49da83d..a8c5d63 100644 --- a/cfspopcon/named_options.py +++ b/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() diff --git a/cfspopcon/variables.yaml b/cfspopcon/variables.yaml index 4a4bf7c..5abd505 100644 --- a/cfspopcon/variables.yaml +++ b/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: diff --git a/docs/conf.py b/docs/conf.py index 5377c57..71aa91d 100644 --- a/docs/conf.py +++ b/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 diff --git a/tests/test_scrape_off_layer_model.py b/tests/test_scrape_off_layer_model.py new file mode 100644 index 0000000..d256218 --- /dev/null +++ b/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)