Implement model to calculate edge impurity conc.

cfspopcon/algorithms/__init__.py

@@ -6,6 +6,7 @@
 from .beta import calc_beta
 from .composite_algorithm import predictive_popcon
 from .core_radiated_power import calc_core_radiated_power
+from .edge_impurity_concentration import calc_edge_impurity_concentration
 from .extrinsic_core_radiator import calc_extrinsic_core_radiator
 from .fusion_gain import calc_fusion_gain
 from .geometry import calc_geometry
@@ -40,6 +41,7 @@
     Algorithms["two_point_model_fixed_tet"]: two_point_model_fixed_tet,
     Algorithms["calc_zeff_and_dilution_from_impurities"]: calc_zeff_and_dilution_from_impurities,
     Algorithms["use_LOC_tau_e_below_threshold"]: use_LOC_tau_e_below_threshold,
+    Algorithms["calc_edge_impurity_concentration"]: calc_edge_impurity_concentration,
     **SINGLE_FUNCTIONS,
 }
 
@@ -53,22 +55,6 @@ def get_algorithm(algorithm: Union[Algorithms, str]) -> Union[Algorithm, Composi
 
 
 __all__ = [
-    "calc_beta",
-    "calc_core_radiated_power",
-    "calc_extrinsic_core_radiator",
-    "calc_fusion_gain",
-    "calc_geometry",
-    "calc_heat_exhaust",
-    "calc_ohmic_power",
-    "calc_peaked_profiles",
-    "calc_plasma_current_from_q_star",
-    "calc_power_balance_from_tau_e",
-    "predictive_popcon",
-    "calc_q_star_from_plasma_current",
-    "two_point_model_fixed_fpow",
-    "two_point_model_fixed_qpart",
-    "two_point_model_fixed_tet",
-    "calc_zeff_and_dilution_from_impurities",
     "ALGORITHMS",
     "get_algorithm",
 ]

cfspopcon/algorithms/edge_impurity_concentration.py

@@ -0,0 +1,71 @@
+"""Run the two point model with a fixed sheath entrance temperature."""
+
+
+from ..atomic_data import read_atomic_data
+from ..formulas.scrape_off_layer_model import build_L_int_integrator, calc_required_edge_impurity_concentration
+from ..named_options import Impurity
+from ..unit_handling import Unitfull, convert_to_default_units, ureg
+from .algorithm_class import Algorithm
+
+RETURN_KEYS = [
+    "edge_impurity_concentration",
+]
+
+
+def run_calc_edge_impurity_concentration(
+    edge_impurity_species: Impurity,
+    q_parallel: Unitfull,
+    SOL_power_loss_fraction: Unitfull,
+    target_electron_temp: Unitfull,
+    upstream_electron_temp: Unitfull,
+    upstream_electron_density: Unitfull,
+    kappa_e0: Unitfull,
+    lengyel_overestimation_factor: Unitfull,
+    reference_electron_density: Unitfull = 1.0 * ureg.n20,
+    reference_ne_tau: Unitfull = 1.0 * ureg.n20 * ureg.ms,
+) -> dict[str, Unitfull]:
+    """Calculate the impurity concentration required to cool the scrape-off-layer using the Lengyel model.
+
+    Args:
+        edge_impurity_species: :term:`glossary link<edge_impurity_species>`
+        reference_electron_density: :term:`glossary link<reference_electron_density>`
+        reference_ne_tau: :term:`glossary link<reference_ne_tau>`
+        q_parallel: :term:`glossary link<q_parallel>`
+        SOL_power_loss_fraction: :term:`glossary link<SOL_power_loss_fraction>`
+        target_electron_temp: :term:`glossary link<target_electron_temp>`
+        upstream_electron_temp: :term:`glossary link<upstream_electron_temp>`
+        upstream_electron_density: :term:`glossary link<upstream_electron_density>`
+        kappa_e0: :term:`glossary link<kappa_e0>`
+        lengyel_overestimation_factor: :term:`glossary link<lengyel_overestimation_factor>`
+
+    Returns:
+        :term:`edge_impurity_concentration`
+    """
+    atomic_data = read_atomic_data()
+
+    L_int_integrator = build_L_int_integrator(
+        atomic_data=atomic_data,
+        impurity_species=edge_impurity_species,
+        reference_electron_density=reference_electron_density,
+        reference_ne_tau=reference_ne_tau,
+    )
+
+    edge_impurity_concentration = calc_required_edge_impurity_concentration(
+        L_int_integrator=L_int_integrator,
+        q_parallel=q_parallel,
+        SOL_power_loss_fraction=SOL_power_loss_fraction,
+        target_electron_temp=target_electron_temp,
+        upstream_electron_temp=upstream_electron_temp,
+        upstream_electron_density=upstream_electron_density,
+        kappa_e0=kappa_e0,
+        lengyel_overestimation_factor=lengyel_overestimation_factor,
+    )
+
+    local_vars = locals()
+    return {key: convert_to_default_units(local_vars[key], key) for key in RETURN_KEYS}
+
+
+calc_edge_impurity_concentration = Algorithm(
+    function=run_calc_edge_impurity_concentration,
+    return_keys=RETURN_KEYS,
+)

cfspopcon/algorithms/single_functions.py

@@ -71,5 +71,10 @@
     return_keys=["plasma_stored_energy"],
     name="calc_plasma_stored_energy",
 )
+calc_upstream_electron_density = Algorithm.from_single_function(
+    lambda nesep_over_nebar, average_electron_density: nesep_over_nebar * average_electron_density,
+    return_keys=["upstream_electron_density"],
+    name="calc_upstream_electron_density",
+)
 
 SINGLE_FUNCTIONS = {Algorithms[key]: val for key, val in locals().items() if isinstance(val, Algorithm)}

cfspopcon/algorithms/two_point_model_fixed_tet.py

@@ -21,8 +21,7 @@ def run_two_point_model_fixed_tet(
     target_electron_temp: Unitfull,
     q_parallel: Unitfull,
     parallel_connection_length: Unitfull,
-    average_electron_density: Unitfull,
-    nesep_over_nebar: Unitfull,
+    upstream_electron_density: Unitfull,
     toroidal_flux_expansion: Unitfull,
     fuel_average_mass_number: Unitfull,
     kappa_e0: Unitfull,
@@ -35,7 +34,7 @@ def run_two_point_model_fixed_tet(
         q_parallel: :term:`glossary link<q_parallel>`
         parallel_connection_length: :term:`glossary link<parallel_connection_length>`
         average_electron_density: :term:`glossary link<average_electron_density>`
-        nesep_over_nebar: :term:`glossary link<nesep_over_nebar>`
+        upstream_electron_density: :term:`glossary link<upstream_electron_density>`
         toroidal_flux_expansion: :term:`glossary link<toroidal_flux_expansion>`
         fuel_average_mass_number: :term:`glossary link<fuel_average_mass_number>`
         kappa_e0: :term:`glossary link<kappa_e0>`
@@ -48,7 +47,7 @@ def run_two_point_model_fixed_tet(
         target_electron_temp=target_electron_temp,
         parallel_heat_flux_density=q_parallel,
         parallel_connection_length=parallel_connection_length,
-        upstream_electron_density=nesep_over_nebar * average_electron_density,
+        upstream_electron_density=upstream_electron_density,
         toroidal_flux_expansion=toroidal_flux_expansion,
         fuel_average_mass_number=fuel_average_mass_number,
         kappa_e0=kappa_e0,

cfspopcon/formulas/__init__.py

@@ -1,6 +1,6 @@
 """Formulas for POPCONs analysis."""
 
-from . import energy_confinement_time_scalings, fusion_reaction_data, plasma_profile_data, radiated_power
+from . import energy_confinement_time_scalings, fusion_reaction_data, plasma_profile_data, radiated_power, scrape_off_layer_model
 from .average_fuel_ion_mass import calc_fuel_average_mass_number
 from .beta import calc_beta_normalised, calc_beta_poloidal, calc_beta_toroidal, calc_beta_total
 from .confinement_regime_threshold_powers import (
@@ -41,54 +41,55 @@
 )
 
 __all__ = [
-    "energy_confinement_time_scalings",
-    "fusion_reaction_data",
-    "calc_density_peaking",
-    "calc_effective_collisionality",
-    "plasma_profile_data",
-    "calc_fuel_average_mass_number",
+    "calc_1D_plasma_profiles",
     "calc_B_pol_omp",
     "calc_B_tor_omp",
     "calc_beta_normalised",
     "calc_beta_poloidal",
     "calc_beta_toroidal",
-    "calc_troyon_limit",
-    "calc_greenwald_density_limit",
     "calc_beta_total",
     "calc_bootstrap_fraction",
+    "calc_bremsstrahlung_radiation",
+    "calc_change_in_dilution",
+    "calc_change_in_zeff",
+    "calc_confinement_transition_threshold_power",
     "calc_current_relaxation_time",
+    "calc_density_peaking",
+    "calc_effective_collisionality",
     "calc_f_shaping",
+    "calc_fuel_average_mass_number",
     "calc_fusion_power",
+    "calc_greenwald_density_limit",
     "calc_greenwald_fraction",
+    "calc_impurity_charge_state",
+    "calc_impurity_radiated_power_mavrin_coronal",
+    "calc_impurity_radiated_power_mavrin_noncoronal",
+    "calc_impurity_radiated_power_post_and_jensen",
+    "calc_impurity_radiated_power_radas",
+    "calc_impurity_radiated_power",
     "calc_LH_transition_threshold_power",
     "calc_LI_transition_threshold_power",
     "calc_loop_voltage",
-    "calc_resistivity_trapped_enhancement",
     "calc_neoclassical_loop_resistivity",
     "calc_neutron_flux_to_walls",
     "calc_normalised_collisionality",
-    "calc_1D_plasma_profiles",
     "calc_ohmic_power",
     "calc_peak_pressure",
     "calc_plasma_current",
     "calc_plasma_surface_area",
     "calc_plasma_volume",
+    "calc_q_star",
+    "calc_resistivity_trapped_enhancement",
     "calc_rho_star",
     "calc_Spitzer_loop_resistivity",
-    "calc_q_star",
-    "calc_triple_product",
-    "calc_tau_e_and_P_in_from_scaling",
-    "thermal_calc_gain_factor",
-    "calc_confinement_transition_threshold_power",
-    "calc_change_in_zeff",
-    "calc_change_in_dilution",
-    "calc_bremsstrahlung_radiation",
     "calc_synchrotron_radiation",
-    "calc_impurity_radiated_power_post_and_jensen",
-    "calc_impurity_radiated_power_radas",
-    "calc_impurity_radiated_power",
-    "calc_impurity_radiated_power_mavrin_coronal",
-    "calc_impurity_radiated_power_mavrin_noncoronal",
+    "calc_tau_e_and_P_in_from_scaling",
+    "calc_triple_product",
+    "calc_troyon_limit",
+    "energy_confinement_time_scalings",
+    "fusion_reaction_data",
+    "plasma_profile_data",
     "radiated_power",
-    "calc_impurity_charge_state",
+    "scrape_off_layer_model",
+    "thermal_calc_gain_factor",
 ]

cfspopcon/formulas/scrape_off_layer_model/__init__.py

@@ -2,6 +2,7 @@
 
 These are mostly based on the two-point-model, from :cite:`stangeby_2018`.
 """
+from .edge_impurity_concentration import build_L_int_integrator, calc_required_edge_impurity_concentration
 from .lambda_q import calc_lambda_q
 from .parallel_heat_flux_density import calc_parallel_heat_flux_density
 from .solve_target_first_two_point_model import solve_target_first_two_point_model
@@ -12,4 +13,6 @@
     "calc_parallel_heat_flux_density",
     "calc_lambda_q",
     "solve_target_first_two_point_model",
+    "calc_required_edge_impurity_concentration",
+    "build_L_int_integrator",
 ]

cfspopcon/formulas/scrape_off_layer_model/edge_impurity_concentration.py

@@ -0,0 +1,87 @@
+"""Lengyel model to compute the edge impurity concentration."""
+from typing import Callable
+
+import numpy as np
+import xarray as xr
+from scipy.interpolate import InterpolatedUnivariateSpline
+
+from ...named_options import Impurity
+from ...unit_handling import Unitfull, convert_units, magnitude, ureg, wraps_ufunc
+
+
+def build_L_int_integrator(
+    atomic_data: dict[Impurity, Unitfull],
+    impurity_species: Impurity,
+    reference_electron_density: Unitfull,
+    reference_ne_tau: Unitfull,
+) -> Callable[[Unitfull, Unitfull], Unitfull]:
+    r"""Build an interpolator to calculate L_{int}$ between arbitrary temperature points.
+
+    $L_int = \int_a^b L_z(T_e) sqrt(T_e) dT_e$ where $L_z$ is a cooling curve for an impurity species.
+    This is used in the calculation of the radiated power associated with a given impurity.
+
+    Args:
+        atomic_data: :term:`glossary link<atomic_data>`
+        impurity_species: [] :term:`glossary link<impurity_species>`
+        reference_electron_density: [n20] :term:`glossary link<reference_electron_density>`
+        reference_ne_tau: [n20 * ms] :term:`glossary link<reference_ne_tau>`
+
+    Returns:
+        :term:`L_int_interpolator`
+    """
+    if isinstance(impurity_species, xr.DataArray):
+        impurity_species = impurity_species.item()
+
+    Lz_curve = atomic_data[impurity_species].noncoronal_electron_emission_prefactor.sel(
+        dim_log_electron_density=np.log10(magnitude(convert_units(reference_electron_density, ureg.m**-3))),
+        dim_log_ne_tau=np.log10(magnitude(convert_units(reference_ne_tau, ureg.m**-3 * ureg.s))),
+    )
+
+    log_electron_temp = Lz_curve.dim_log_electron_temperature
+    electron_temp = np.power(10, log_electron_temp)
+    Lz_sqrt_Te = Lz_curve * np.sqrt(electron_temp)
+
+    interpolator = InterpolatedUnivariateSpline(electron_temp, magnitude(Lz_sqrt_Te))
+    def L_int(start_temp, stop_temp):
+        return interpolator.integral(start_temp, stop_temp)
+
+    return wraps_ufunc(
+        input_units=dict(start_temp=ureg.eV, stop_temp=ureg.eV), return_units=dict(L_int=ureg.W * ureg.m**3 * ureg.eV**1.5)
+    )(L_int)
+
+
+def calc_required_edge_impurity_concentration(
+    L_int_integrator: Callable[[Unitfull, Unitfull], Unitfull],
+    q_parallel: Unitfull,
+    SOL_power_loss_fraction: Unitfull,
+    target_electron_temp: Unitfull,
+    upstream_electron_temp: Unitfull,
+    upstream_electron_density: Unitfull,
+    kappa_e0: Unitfull,
+    lengyel_overestimation_factor: Unitfull,
+) -> Unitfull:
+    """Calculate the relative concentration of an edge impurity required to achieve a given SOL power loss fraction.
+
+    N.b. this function does not ensure consistency of the calculated impurity concentration
+    with the parallel temperature profile. You may wish to implement an iterative solver
+    to find a consistent set of L_parallel, T_t and T_u.
+
+    Args:
+        L_int_integrator: :term:`glossary link<L_int_integrator>`
+        q_parallel: :term:`glossary link<q_parallel>`
+        SOL_power_loss_fraction: :term:`glossary link<SOL_power_loss_fraction>`
+        target_electron_temp: :term:`glossary link<target_electron_temp>`
+        upstream_electron_temp: :term:`glossary link<upstream_electron_temp>`
+        upstream_electron_density: :term:`glossary link<upstream_electron_density>`
+        kappa_e0: :term:`glossary link<kappa_e0>`
+        lengyel_overestimation_factor: :term:`glossary link<lengyel_overestimation_factor>`
+
+    Returns:
+        :term:`impurity_concentration`
+    """
+    L_int = L_int_integrator(target_electron_temp, upstream_electron_temp)
+
+    numerator = q_parallel**2 - ((1.0 - SOL_power_loss_fraction) * q_parallel) ** 2
+    denominator = 2.0 * kappa_e0 * (upstream_electron_density * upstream_electron_temp) ** 2 * L_int
+
+    return numerator / denominator / lengyel_overestimation_factor

cfspopcon/helpers.py

@@ -33,6 +33,8 @@ def convert_named_options(key: str, val: Any) -> Any:  # noqa: PLR0911, PLR0912
         return make_impurities_array(list(val.keys()), list(val.values()))
     elif key == "core_radiator":
         return Impurity[val]
+    elif key == "edge_impurity_species":
+        return Impurity[val]
     elif key == "lambda_q_scaling":
         return LambdaQScaling[val]
     elif key == "SOL_momentum_loss_function":

cfspopcon/named_options.py

@@ -40,6 +40,8 @@ class Algorithms(Enum):
     calc_P_SOL = auto()
     use_LOC_tau_e_below_threshold = auto()
     calc_plasma_stored_energy = auto()
+    calc_edge_impurity_concentration = auto()
+    calc_upstream_electron_density = auto()
 
 
 class ProfileForm(Enum):

cfspopcon/unit_handling/default_units.py

@@ -131,6 +131,10 @@
     vertical_minor_radius="m",
     z_effective="",
     zeff_change_from_core_rad="",
+    edge_impurity_species=None,
+    lengyel_overestimation_factor="",
+    upstream_electron_density="n19",
+    edge_impurity_concentration="",
 )
 
 
@@ -252,7 +256,7 @@ def convert_to_default_units(value: Quantity, key: str) -> Quantity:
 def convert_to_default_units(value: Union[float, Quantity, xr.DataArray], key: str) -> Union[float, Quantity, xr.DataArray]:
     """Convert an array or scalar to default units."""
     unit = DEFAULT_UNITS[key]
-    if unit is None or unit == "":
+    if unit is None:
         return value
     elif isinstance(value, (xr.DataArray, Quantity)):
         return convert_units(value, unit)

example_cases/SPARC_PRD/input.yaml

@@ -42,7 +42,9 @@ algorithms:
   - calc_P_SOL
   - calc_average_total_pressure
   - calc_heat_exhaust
+  - calc_upstream_electron_density
   - two_point_model_fixed_tet
+  - calc_edge_impurity_concentration
   # Finally, we calculate several parameters which aren't used in other
   # calculations, but which are useful for characterizing operational
   # points. These can be used later when masking inaccessible operational
@@ -166,3 +168,8 @@ fraction_of_P_SOL_to_divertor: 0.6
 kappa_e0: 2600.0
 # Calculate P_rad_SOL such that the target electron temperature in eV reaches this value (for FixedTargetElectronTemp)
 target_electron_temp: 25.0
+
+# Which species should be injected to cool the scrape-off-layer?
+edge_impurity_species: Argon
+# Factor applied to the result of the Lengyel model to give an absolute low-Z impurity concentration.
+lengyel_overestimation_factor: 4.3