Add option for density pedestal setting

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

@@ -148,4 +148,32 @@ $$\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}
+$$
+
+
+$\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
+
+------
+
 [^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

@@ -60,6 +60,7 @@
 from process.data_structure.tfcoil_variables import init_tfcoil_variables
 from process.data_structure.times_variables import init_times_variables
 from process.data_structure.vacuum_variables import init_vacuum_variables
+from process.models.physics.profiles import DensityProfilePedestalType
 from process.models.stellarator.initialization import st_init
 from process.models.superconductors import (
     SuperconductorMaterial,
@@ -503,10 +504,10 @@ def check_process(inputs):  # noqa: ARG001
         # Density checks
         # Case where pedestal density is set manually
         if (
-            data_structure.physics_variables.f_nd_plasma_pedestal_greenwald < 0
-            or not (
-                data_structure.numerics.ixc[: data_structure.numerics.nvar] == 145
-            ).any()
+            DensityProfilePedestalType(
+                data_structure.physics_variables.i_nd_plasma_pedestal_separatrix
+            )
+            == DensityProfilePedestalType.USER_INPUT
         ):
             # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron but it is less than nd_plasma_separatrix_electron.
             if (

process/core/input.py

@@ -171,10 +171,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.01, 0.9)
     ),
     "f_plasma_fuel_helium3": InputVariable(
         data_structure.physics_variables, float, range=(-1.0, 5.0)
@@ -1083,6 +1083,9 @@ def __post_init__(self):
     "nd_plasma_separatrix_electron": InputVariable(
         data_structure.physics_variables, float, range=(0.0, 1e21)
     ),
+    "i_nd_plasma_pedestal_electron": InputVariable(
+        data_structure.physics_variables, int, choices=[0, 1]
+    ),
     "nflutfmax": InputVariable(
         data_structure.constraint_variables, float, range=(0.0, 1e24)
     ),

process/core/io/plot/summary.py

@@ -4015,7 +4015,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}"
@@ -4261,33 +4264,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/solver/iteration_variables.py

@@ -260,10 +260,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", data_structure.ife_variables, 5.0e2, 3.0e3),
     156: IterationVariable("rrin", data_structure.ife_variables, 1.0, 1.0e1),

process/data_structure/physics_variables.py

@@ -422,16 +422,12 @@
 
 
 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`)
+"""fraction of Greenwald density to set as pedestal-top 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`)
+"""fraction of Greenwald density to set as separatrix density.
 """
 
 
@@ -601,6 +597,12 @@
 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"""
@@ -1207,6 +1209,8 @@
 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`"""
@@ -1488,6 +1492,7 @@ def init_physics_variables():
         ffwal, \
         f_nd_plasma_pedestal_greenwald, \
         f_nd_plasma_separatrix_greenwald, \
+        i_nd_plasma_pedestal_separatrix, \
         f_plasma_fuel_helium3, \
         figmer, \
         fkzohm, \
@@ -1653,6 +1658,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, \
@@ -1777,6 +1783,7 @@ def init_physics_variables():
     ffwal = 0.92
     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
@@ -1944,6 +1951,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

process/models/physics/physics.py

@@ -35,7 +35,10 @@
 from process.models.physics.density_limit import PlasmaDensityLimit
 from process.models.physics.exhaust import PlasmaExhaust
 from process.models.physics.l_h_transition import PlasmaConfinementTransition
-from process.models.physics.profiles import PlasmaProfileShapeType
+from process.models.physics.profiles import (
+    DensityProfilePedestalType,
+    PlasmaProfileShapeType,
+)
 
 if TYPE_CHECKING:
     from process.models.physics.plasma_current import (
@@ -340,24 +343,8 @@ def run(self):
         if (
             PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
             == PlasmaProfileShapeType.PEDESTAL_PROFILE
-        ) and (physics_variables.f_nd_plasma_pedestal_greenwald >= 0e0):
-            physics_variables.nd_plasma_pedestal_electron = (
-                physics_variables.f_nd_plasma_pedestal_greenwald
-                * 1.0e14
-                * physics_variables.plasma_current
-                / (np.pi * physics_variables.rminor * physics_variables.rminor)
-            )
-
-        if (
-            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
-            == PlasmaProfileShapeType.PEDESTAL_PROFILE
-        ) and (physics_variables.f_nd_plasma_separatrix_greenwald >= 0e0):
-            physics_variables.nd_plasma_separatrix_electron = (
-                physics_variables.f_nd_plasma_separatrix_greenwald
-                * 1.0e14
-                * physics_variables.plasma_current
-                / (np.pi * physics_variables.rminor * physics_variables.rminor)
-            )
+        ):
+            self.plasma_profile.neprofile.set_pedestal_and_separatrix_values()
 
         self.plasma_profile.run()
 
@@ -2069,6 +2056,12 @@ def outplas(self):
             "(temp_plasma_electron_vol_avg_kev)",
             physics_variables.temp_plasma_electron_vol_avg_kev,
         )
+        po.ovarrf(
+            self.outfile,
+            "Line averaged electron temperature (keV)",
+            "(temp_plasma_electron_line_avg_kev)",
+            physics_variables.temp_plasma_electron_line_avg_kev,
+        )
         po.ovarrf(
             self.outfile,
             "Ratio of ion to electron volume-averaged temperature",
@@ -2414,7 +2407,12 @@ def outplas(self):
                 "(radius_plasma_pedestal_density_norm)",
                 physics_variables.radius_plasma_pedestal_density_norm,
             )
-            if physics_variables.f_nd_plasma_pedestal_greenwald >= 0e0:
+            if (
+                DensityProfilePedestalType(
+                    physics_variables.i_nd_plasma_pedestal_separatrix
+                )
+                == DensityProfilePedestalType.GREENWALD_FRACTION
+            ):
                 po.ovarre(
                     self.outfile,
                     "Electron density pedestal height (/m3)",

process/models/physics/plasma_profiles.py

@@ -127,6 +127,14 @@ def parabolic_paramterisation(self):
             / sp.special.gamma(physics_variables.alphan + 1.5)
         )
 
+        physics_variables.temp_plasma_electron_line_avg_kev = (
+            physics_variables.temp_plasma_electron_vol_avg_kev
+            * (1.0 + physics_variables.alphat)
+            * (sp.special.gamma(0.5) / 2.0)
+            * sp.special.gamma(physics_variables.alphat + 1.0)
+            / sp.special.gamma(physics_variables.alphat + 1.5)
+        )
+
         #  Density-weighted temperatures
 
         physics_variables.temp_plasma_electron_density_weighted_kev = (
@@ -203,11 +211,15 @@ def pedestal_parameterisation(self):
             / physics_variables.temp_plasma_electron_vol_avg_kev
         )
 
-        #  Line-averaged electron density
+        #  Line-averaged electron density and temperature
         #  = integral(n(rho).drho)
 
         physics_variables.nd_plasma_electron_line = self.neprofile.profile_integ
 
+        physics_variables.temp_plasma_electron_line_avg_kev = (
+            self.teprofile.profile_integ
+        )
+
         #  Scrape-off density / volume averaged density
         #  (Input value is used if i_plasma_pedestal = 0)
 

process/models/physics/profiles.py

@@ -6,6 +6,7 @@
 import scipy as sp
 
 from process.data_structure import physics_variables
+from process.models.physics.density_limit import PlasmaDensityLimit
 
 logger = logging.getLogger(__name__)
 
@@ -83,6 +84,13 @@ def integrate_profile_y(self):
         )
 
 
+class DensityProfilePedestalType(IntEnum):
+    """Enum for i_nd_plasma_pedestal_separatrix types"""
+
+    USER_INPUT = 0
+    GREENWALD_FRACTION = 1
+
+
 class NeProfile(Profile):
     """Electron density profile class. Contains a function to calculate the electron density profile and
     store the data.
@@ -220,6 +228,33 @@ def ncore(
             ncore = 1.0e-6
         return ncore
 
+    def set_pedestal_and_separatrix_values(self):
+        """Sets the pedestal and separatrix density values based on the user input or greenwald fraction method."""
+
+        if (
+            DensityProfilePedestalType(physics_variables.i_nd_plasma_pedestal_separatrix)
+            == DensityProfilePedestalType.USER_INPUT
+        ):
+            pass
+        elif (
+            DensityProfilePedestalType(physics_variables.i_nd_plasma_pedestal_separatrix)
+            == DensityProfilePedestalType.GREENWALD_FRACTION
+        ):
+            physics_variables.nd_plasma_pedestal_electron = (
+                physics_variables.f_nd_plasma_pedestal_greenwald
+                * PlasmaDensityLimit.calculate_greenwald_density_limit(
+                    c_plasma=physics_variables.plasma_current,
+                    rminor=physics_variables.rminor,
+                )
+            )
+            physics_variables.nd_plasma_separatrix_electron = (
+                physics_variables.f_nd_plasma_separatrix_greenwald
+                * PlasmaDensityLimit.calculate_greenwald_density_limit(
+                    c_plasma=physics_variables.plasma_current,
+                    rminor=physics_variables.rminor,
+                )
+            )
+
     def set_physics_variables(self):
         """Calculates and sets physics variables required for the profile."""