Add option for density pedestal setting

documentation/source/physics-models/profiles/plasma_density_profile.md

@@ -148,4 +148,31 @@ $$\begin{aligned}
 
 5. Profile is then integrated with `integrate_profile_y()` using Simpsons integration from the profile abstract base class
 
+-----------------------
+
+### Setting pedestal and separatrix values | `set_pedestal_and_separatrix_values()`
+
+The switch `i_nd_plasma_pedestal_separatrix` controls how the values of the density pedestal and separatrix are set.
+
+#### User input
+
+If `i_nd_plasma_pedestal_separatrix == 0` then the values of `nd_plasma_pedestal_electron` and `nd_plasma_separatrix_electron` are taken directly from the input file
+
+#### Fraction of Greenwald Limit
+
+If `i_nd_plasma_pedestal_separatrix == 1`, the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ are set as fractions of the [Greenwald](https://wiki.fusion.ciemat.es/wiki/Greenwald_limit) limit such as:
+
+$$
+n_{\text{ped}} = \overbrace{f_{\text{GW,ped}}}^{\texttt{f_nd_plasma_pedestal_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
+$$
+
+$$
+n_{\text{sep}} = \overbrace{f_{\text{GW,sep}}}^{\texttt{f_nd_plasma_separatrix_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
+$$
+
+
+`f_nd_plasma_pedestal_greenwald` and `f_nd_plasma_separatrix_greenwald` can be set as iteration variables respectively by using `ixc = 145`
+and `ixc = 152` respectively
+
+
 [^1]: Jean, J. (2011). *HELIOS: A Zero-Dimensional Tool for Next Step and Reactor Studies*. Fusion Science and Technology, 59(2), 308–349. <https://doi.org/10.13182/FST11-A11650>

documentation/source/physics-models/profiles/plasma_profiles.md

@@ -636,24 +636,6 @@ The same function is run from the `i_plasma_pedestal == 0 ` profile case, found
 
 --------
 
-### Setting pedestal values as fractions of the Greenwald limit
-
-By default, the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ are set as fractions of the [Greenwald](https://wiki.fusion.ciemat.es/wiki/Greenwald_limit) limit such as:
-
-$$
-n_{\text{ped}} = \overbrace{f_{\text{GW,ped}}}^{\texttt{f_nd_plasma_pedestal_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
-$$
-
-$$
-n_{\text{sep}} = \overbrace{f_{\text{GW,sep}}}^{\texttt{f_nd_plasma_separatrix_greenwald}} \times \frac{I_p [\text{A}]}{\pi a^2 [\text{m}^2]} \times 10^{14}
-$$
-
-To set the values of $n_{\text{ped}}$ and $n_{\text{sep}}$ directly, the user can input the value of $\texttt{f_nd_plasma_pedestal_greenwald}$ or $\texttt{f_nd_plasma_separatrix_greenwald}$ to be less than 0.0 (i.e negative) to prevent the Greenwald fraction value being set.
-
-$\texttt{f_nd_plasma_pedestal_greenwald}$ and $\texttt{f_nd_plasma_separatrix_greenwald}$ can be set as iteration variables respectively by using `ixc = 45`
-and `ixc = 152` respectively
-
-------
 
 ### Pedestal Density Upper limit
 

process/core/init.py

@@ -30,6 +30,7 @@
     init_superconducting_tf_coil_variables,
 )
 from process.data_structure.tfcoil_variables import init_tfcoil_variables
+from process.models.physics.profiles import DensityProfilePedestalType
 from process.models.stellarator.initialization import st_init
 from process.models.superconductors import (
     SuperconductorMaterial,
@@ -443,45 +444,53 @@ def check_process(inputs, data):  # noqa: ARG001
             )
 
         # Density checks
-        # Case where pedestal density is set manually
+        # Issue #589: Pedestal density is lower than separatrix density
+        pedestal_type = DensityProfilePedestalType(
+            data_structure.physics_variables.i_nd_plasma_pedestal_separatrix
+        )
         if (
-            data_structure.physics_variables.f_nd_plasma_pedestal_greenwald < 0
-            or not (
-                data_structure.numerics.ixc[: data_structure.numerics.nvar] == 145
-            ).any()
+            pedestal_type == DensityProfilePedestalType.USER_INPUT
+            and data_structure.physics_variables.nd_plasma_pedestal_electron
+            < data_structure.physics_variables.nd_plasma_separatrix_electron
+        ) or (
+            pedestal_type == DensityProfilePedestalType.GREENWALD_FRACTION
+            and data_structure.physics_variables.f_nd_plasma_pedestal_greenwald
+            < data_structure.physics_variables.f_nd_plasma_separatrix_greenwald
         ):
-            # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron but it is less than nd_plasma_separatrix_electron.
-            if (
-                data_structure.physics_variables.nd_plasma_pedestal_electron
-                < data_structure.physics_variables.nd_plasma_separatrix_electron
-            ):
-                raise ProcessValidationError(
-                    "Density pedestal is lower than separatrix density",
-                    nd_plasma_pedestal_electron=data_structure.physics_variables.nd_plasma_pedestal_electron,
-                    nd_plasma_separatrix_electron=data_structure.physics_variables.nd_plasma_separatrix_electron,
-                )
+            raise ProcessValidationError(
+                "Density pedestal is lower than separatrix density",
+                **(
+                    {
+                        "nd_plasma_pedestal_electron": data_structure.physics_variables.nd_plasma_pedestal_electron,
+                        "nd_plasma_separatrix_electron": data_structure.physics_variables.nd_plasma_separatrix_electron,
+                    }
+                    if pedestal_type == DensityProfilePedestalType.USER_INPUT
+                    else {
+                        "f_nd_plasma_pedestal_greenwald": data_structure.physics_variables.f_nd_plasma_pedestal_greenwald,
+                        "f_nd_plasma_separatrix_greenwald": data_structure.physics_variables.f_nd_plasma_separatrix_greenwald,
+                    }
+                ),
+            )
 
-            # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron,
-            # but pedestal width = 0.
-            if (
-                abs(
-                    data_structure.physics_variables.radius_plasma_pedestal_density_norm
-                    - 1.0
-                )
-                <= 1e-7
-                and (
-                    data_structure.physics_variables.nd_plasma_pedestal_electron
-                    - data_structure.physics_variables.nd_plasma_separatrix_electron
-                )
-                >= 1e-7
-            ):
-                warn(
-                    "Density pedestal is at plasma edge "
-                    f"({data_structure.physics_variables.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron "
-                    f"({data_structure.physics_variables.nd_plasma_pedestal_electron}) differs from "
-                    f"nd_plasma_separatrix_electron ({data_structure.physics_variables.nd_plasma_separatrix_electron})",
-                    stacklevel=2,
-                )
+        if (
+            abs(
+                data_structure.physics_variables.radius_plasma_pedestal_density_norm
+                - 1.0
+            )
+            <= 1e-7
+            and (
+                data_structure.physics_variables.nd_plasma_pedestal_electron
+                - data_structure.physics_variables.nd_plasma_separatrix_electron
+            )
+            >= 1e-7
+        ):
+            warn(
+                "Density pedestal is at plasma edge "
+                f"({data_structure.physics_variables.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron "
+                f"({data_structure.physics_variables.nd_plasma_pedestal_electron}) differs from "
+                f"nd_plasma_separatrix_electron ({data_structure.physics_variables.nd_plasma_separatrix_electron})",
+                stacklevel=2,
+            )
 
         # Issue #862 : Variable nd_plasma_electron_on_axis/nd_plasma_pedestal_electron ratio without constraint eq 81 (nd_plasma_electron_on_axis>nd_plasma_pedestal_electron)
         #  -> Potential hollowed density profile

process/core/input.py

@@ -176,10 +176,10 @@ def __post_init__(self):
     ),
     "ffwal": InputVariable(data_structure.physics_variables, float, range=(0.0, 10.0)),
     "f_nd_plasma_pedestal_greenwald": InputVariable(
-        data_structure.physics_variables, float, range=(-1.0, 5.0)
+        data_structure.physics_variables, float, range=(0.1, 1.5)
     ),
     "f_nd_plasma_separatrix_greenwald": InputVariable(
-        data_structure.physics_variables, float, range=(-1.0, 5.0)
+        data_structure.physics_variables, float, range=(0.001, 0.9)
     ),
     "f_plasma_fuel_helium3": InputVariable(
         data_structure.physics_variables, float, range=(-1.0, 5.0)
@@ -692,6 +692,9 @@ def __post_init__(self):
     "nd_plasma_separatrix_electron": InputVariable(
         data_structure.physics_variables, float, range=(0.0, 1e21)
     ),
+    "i_nd_plasma_pedestal_separatrix": InputVariable(
+        data_structure.physics_variables, int, choices=[0, 1]
+    ),
     "nflutfmax": InputVariable("constraints", float, range=(0.0, 1e24)),
     "oacdcp": InputVariable(
         data_structure.tfcoil_variables, float, range=(10000.0, 1000000000.0)

process/core/io/mfile/comparison.py

@@ -117,7 +117,7 @@
     "temp_plasma_electron_on_axis_kev",
     "nd_plasma_electrons_vol_avg",
     "nd_plasma_electron_on_axis",
-    "dnla_gw",
+    "f_nd_plasma_greenwald",
     "temp_plasma_separatrix_kev",
     "nd_plasma_separatrix_electron",
     "temp_plasma_pedestal_kev",

process/core/io/plot/summary.py

@@ -3163,7 +3163,7 @@ def plot_main_plasma_information(
         f"$\\mathbf{{Density \\ limit:}}$\n"
         f"({DensityLimitModel(int(mfile.get('i_density_limit', scan=scan))).full_name})\n"
         f"$n_{{\\text{{e,limit}}}}: {mfile.get('nd_plasma_electrons_max', scan=scan):.3e} \\ m^{{-3}}$\n"
-        f"$f_{{\\text{{GW}}}}$: {mfile.get('dnla_gw', scan=scan):.4f}"
+        f"$f_{{\\text{{GW}}}}$: {mfile.get('f_nd_plasma_greenwald', scan=scan):.4f}"
     )
 
     axis.text(
@@ -3857,8 +3857,12 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
     nd_ions_total = mfile.get("nd_plasma_ions_total_vol_avg", scan=scan)
     nd_fuel_ions = mfile.get("nd_plasma_fuel_ions_vol_avg", scan=scan)
     alphan = mfile.get("alphan", scan=scan)
-    fgwped_out = mfile.get("fgwped_out", scan=scan)
-    fgwsep_out = mfile.get("fgwsep_out", scan=scan)
+    f_nd_plasma_pedestal_greenwald = mfile.get(
+        "f_nd_plasma_pedestal_greenwald", scan=scan
+    )
+    f_nd_plasma_separatrix_greenwald = mfile.get(
+        "f_nd_plasma_separatrix_greenwald", scan=scan
+    )
     nd_plasma_electrons_vol_avg = mfile.get("nd_plasma_electrons_vol_avg", scan=scan)
     # find impurity densities
     imp_frac = np.array([
@@ -4048,7 +4052,10 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
 
     # Add text box with density profile parameters
     textstr_density = "\n".join((
-        rf"$\langle n_{{\text{{e}}}} \rangle$: {nd_plasma_electrons_vol_avg:.3e} m$^{{-3}}$",
+        (
+            rf"$\langle n_{{\text{{e}}}} \rangle$: {nd_plasma_electrons_vol_avg:.3e} m$^{{-3}}$"
+            rf"$\hspace{{4}} \overline{{n_{{e}}}}$: {mfile.get('nd_plasma_electron_line', scan=scan):.3e} m$^{{-3}}$"
+        ),
         (
             rf"$n_{{\text{{e,0}}}}$: {ne0:.3e} m$^{{-3}}$"
             rf"$\hspace{{4}} \alpha_{{\text{{n}}}}$: {alphan:.3f}"
@@ -4059,7 +4066,7 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
             f"{nd_fuel_ions / nd_plasma_electrons_vol_avg:.3f}"
         ),
         (
-            rf"$f_{{\text{{GW e,ped}}}}$: {fgwped_out:.3f}"
+            rf"$f_{{\text{{GW e,ped}}}}$: {f_nd_plasma_pedestal_greenwald:.3f}"
             r"$ \hspace{7} \frac{n_{e,0}}{\langle n_e \rangle}$: "
             f"{ne0 / nd_plasma_electrons_vol_avg:.3f}"
         ),
@@ -4069,12 +4076,12 @@ def plot_n_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
             f"{mfile.get('nd_plasma_electron_line', scan=scan) / mfile.get('nd_plasma_electron_max_array(7)', scan=scan):.3f}"
         ),
         rf"$n_{{\text{{e,sep}}}}$: {nd_plasma_separatrix_electron:.3e} m$^{{-3}}$",
-        rf"$f_{{\text{{GW e,sep}}}}$: {fgwsep_out:.3f}",
+        rf"$f_{{\text{{GW e,sep}}}}$: {f_nd_plasma_separatrix_greenwald:.3f}",
     ))
 
     props_density = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5}
     ax_main.text(
-        0.0,
+        -0.05,
         -0.175,
         textstr_density,
         transform=ax_impurity.transAxes,
@@ -4293,33 +4300,34 @@ def plot_t_profiles(prof, demo_ranges: bool, mfile: MFile, scan: int):
     # Add text box with temperature profile parameters
     textstr_temperature = "\n".join((
         (
-            rf"$\langle T_{{\text{{e}}}} \rangle_\text{{V}}$: {mfile.get('temp_plasma_electron_vol_avg_kev', scan=scan):.3f} keV"
-            rf"$\hspace{{3}} \langle T_{{\text{{e}}}} \rangle_\text{{n}}$: {mfile.get('temp_plasma_electron_density_weighted_kev', scan=scan):.3f} keV"
+            rf"$\langle T_{{\text{{e}}}} \rangle_\text{{V}}$:  {mfile.get('temp_plasma_electron_vol_avg_kev', scan=scan):.3f} keV"
+            rf"$\hspace{{2}} \langle T_{{\text{{e}}}} \rangle_\text{{n}}$: {mfile.get('temp_plasma_electron_density_weighted_kev', scan=scan):.3f} keV"
+            rf"$\hspace{{2}} \overline{{T_{{e}}}}$: {mfile.get('temp_plasma_electron_line_avg_kev', scan=scan):.3f} keV"
         ),
         (
-            rf"$T_{{\text{{e,0}}}}$: {te0:.3f} keV"
-            rf"$\hspace{{4}} \alpha_{{\text{{T}}}}$: {alphat:.3f}"
+            rf"$T_{{\text{{e,0}}}}$:    {te0:.3f} keV"
+            rf"$\hspace{{2}} \alpha_{{\text{{T}}}}$:   {alphat:.3f}"
         ),
         (
             rf"$T_{{\text{{e,ped}}}}$: {temp_plasma_pedestal_kev:.3f} keV"
-            r"$ \hspace{4} \frac{\langle T_i \rangle}{\langle T_e \rangle}$: "
+            r"$ \hspace{3} \frac{\langle T_i \rangle}{\langle T_e \rangle}$: "
             f"{f_temp_plasma_ion_electron:.3f}"
         ),
         (
             rf"$\rho_{{\text{{ped,T}}}}$: {radius_plasma_pedestal_temp_norm:.3f}"
-            r"$ \hspace{6} \frac{T_{e,0}}{\langle T_e \rangle}$: "
+            r"$ \hspace{5} \frac{T_{e,0}}{\langle T_e \rangle}$: "
             f"{te0 / te:.3f}"
         ),
         (
             rf"$T_{{\text{{e,sep}}}}$: {temp_plasma_separatrix_kev:.3f} keV"
-            r"$ \hspace{4} \frac{{{\langle T_e \rangle_n}}}{{{\langle T_e \rangle_V}}}$: "
+            r"$ \hspace{3} \frac{{{\langle T_e \rangle_n}}}{{{\langle T_e \rangle_V}}}$: "
             f"{mfile.get('f_temp_plasma_electron_density_vol_avg', scan=scan):.3f}"
         ),
     ))
 
     props_temperature = {"boxstyle": "round", "facecolor": "wheat", "alpha": 0.5}
     prof.text(
-        0.0,
+        -0.1,
         -0.125,
         textstr_temperature,
         transform=prof.transAxes,

process/core/io/variable_metadata.py

@@ -187,7 +187,7 @@ class VariableMetadata:
         description="Average electron density",
         units="m^{-3}",
     ),
-    "dnla_gw": VariableMetadata(
+    "f_nd_plasma_greenwald": VariableMetadata(
         latex=r"$f_{\mathrm{GW}}$", description="Greenwald fraction", units=""
     ),
     "normalised_toroidal_beta": VariableMetadata(

process/core/solver/iteration_variables.py

@@ -235,10 +235,10 @@ class IterationVariable:
         1.0e20,
     ),
     145: IterationVariable(
-        "f_nd_plasma_pedestal_greenwald", data_structure.physics_variables, 0.1, 0.9
+        "f_nd_plasma_pedestal_greenwald", data_structure.physics_variables, 0.1, 1.5
     ),
     152: IterationVariable(
-        "f_nd_plasma_separatrix_greenwald", data_structure.physics_variables, 0.001, 0.5
+        "f_nd_plasma_separatrix_greenwald", data_structure.physics_variables, 0.001, 0.9
     ),
     155: IterationVariable("pfusife", "ife", 5.0e2, 3.0e3),
     156: IterationVariable("rrin", "ife", 1.0, 1.0e1),

process/data_structure/physics_variables.py

@@ -430,18 +430,17 @@ class DivertorNumberModels(IntEnum):
 load calculation (`i_pflux_fw_neutron=1`)
 """
 
+f_nd_plasma_greenwald: float = None
+"""Greenwald fraction of the line averaged electron density. The classic Greenwald
+limit value"""
 
 f_nd_plasma_pedestal_greenwald: float = None
-"""fraction of Greenwald density to set as pedestal-top density. If `<0`, pedestal-top
-density set manually using nd_plasma_pedestal_electron (`i_plasma_pedestal==1`).
-(`iteration variable 145`)
+"""Greenwald fraction of the pedestal density
 """
 
 
 f_nd_plasma_separatrix_greenwald: float = None
-"""fraction of Greenwald density to set as separatrix density. If `<0`, separatrix
-density set manually using nd_plasma_separatrix_electron (`i_plasma_pedestal==1`).
-(`iteration variable 152`)
+"""Greenwald fraction of the separatrix density
 """
 
 
@@ -611,6 +610,12 @@ class DivertorNumberModels(IntEnum):
 nd_plasma_separatrix_electron: float = None
 """electron density at separatrix [m-3] (`i_plasma_pedestal==1)"""
 
+i_nd_plasma_pedestal_separatrix: int = None
+"""switch for pedestal and separatrix density calculation:
+- =0 User input pedestal and separatrix density
+- =1 Calculate pedestal and separatrix density as fraction of Greenwald limit (see `f_nd_plasma_pedestal_greenwald` and `f_nd_plasma_separatrix_greenwald`)
+"""
+
 
 alpha_crit: float = None
 """critical ballooning parameter value"""
@@ -1217,6 +1222,8 @@ class DivertorNumberModels(IntEnum):
 temp_plasma_electron_density_weighted_kev: float = None
 """density weighted average electron temperature (keV)"""
 
+temp_plasma_electron_line_avg_kev: float = None
+"""Line averaged electron temperature (keV)"""
 
 temp_plasma_ion_vol_avg_kev: float = None
 """volume averaged ion temperature (keV). N.B. calculated from temp_plasma_electron_vol_avg_kev if `f_temp_plasma_ion_electron > 0.0`"""
@@ -1628,8 +1635,10 @@ def init_physics_variables():
         f_plasma_fuel_deuterium, \
         f_p_div_lower, \
         ffwal, \
+        f_nd_plasma_greenwald, \
         f_nd_plasma_pedestal_greenwald, \
         f_nd_plasma_separatrix_greenwald, \
+        i_nd_plasma_pedestal_separatrix, \
         f_plasma_fuel_helium3, \
         figmer, \
         fkzohm, \
@@ -1796,6 +1805,7 @@ def init_physics_variables():
         temp_plasma_electron_vol_avg_kev, \
         temp_plasma_electron_on_axis_kev, \
         temp_plasma_electron_density_weighted_kev, \
+        temp_plasma_electron_line_avg_kev, \
         temp_plasma_ion_vol_avg_kev, \
         temp_plasma_ion_on_axis_kev, \
         temp_plasma_ion_density_weighted_kev, \
@@ -1959,8 +1969,10 @@ def init_physics_variables():
     f_plasma_fuel_deuterium = 0.5
     f_p_div_lower = 1.0
     ffwal = 0.92
+    f_nd_plasma_greenwald = 1.0
     f_nd_plasma_pedestal_greenwald = 0.85
     f_nd_plasma_separatrix_greenwald = 0.50
+    i_nd_plasma_pedestal_separatrix = 1
     f_plasma_fuel_helium3 = 0.0
     figmer = 0.0
     fkzohm = 1.0
@@ -2129,6 +2141,7 @@ def init_physics_variables():
     temp_plasma_electron_vol_avg_kev = 12.9
     temp_plasma_electron_on_axis_kev = 0.0
     temp_plasma_electron_density_weighted_kev = 0.0
+    temp_plasma_electron_line_avg_kev = 0.0
     temp_plasma_ion_vol_avg_kev = 12.9
     temp_plasma_ion_on_axis_kev = 0.0
     temp_plasma_ion_density_weighted_kev = 0.0