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

@@ -277,10 +277,10 @@ $$
 
 --------
 
-The line averaged density is then calculated for the profile paramaters
+The line averaged density/temperature is then calculated for the profile paramaters
 
-###### Line averaged density derivation
-Line averaged electron density is calculated by integrating the profile across the normalised width of the profile and then dividing by the width of the integration bounds
+###### Line averaged density/temperature derivation
+Line averaged electron density/temperature is calculated by integrating the profile across the normalised width of the profile and then dividing by the width of the integration bounds. Using the density as an example:
 
 $$
 \overbrace{\bar{n_{\text{e}}}}^{\texttt{nd_plasma_electron_line}} = \frac{\int^1_0 n_0(1-\rho^2)^{\alpha_n} \ d\rho}{\rho}
@@ -612,12 +612,16 @@ $$
 \texttt{f_temp_plasma_electron_density_vol_avg} =\frac{\langle T_{\text{e}} \rangle_{\text{n}}}{\underbrace{\langle T_{\text{e}} \rangle_{\text{V}}}_{\texttt{temp_plasma_electrons_vol_avg}}}
 $$
 
-Calculate the line averaged electron density by integrating the normalised profile using the class [`integrate_profile_y()`](./plasma_profiles_abstract_class.md#calculate-the-profile-integral-value-integrate_profile_y) function
+Calculate the line averaged electron density and temperature by integrating the normalised profile using the class [`integrate_profile_y()`](./plasma_profiles_abstract_class.md#calculate-the-profile-integral-value-integrate_profile_y) function
 
 $$
  \overbrace{\bar{n_{\text{e}}}}^{\texttt{nd_plasma_electron_line}} = \int_0^1{n(\rho) \ d\rho}
 $$
 
+$$
+ \overbrace{\bar{T_{\text{e}}}}^{\texttt{temp_plasma_electron_line_avg_kev}} = \int_0^1{T(\rho) \ d\rho}
+$$
+
 A divertor variable `prn1` is set to be equal to the separatrix density over the mean density:
 
 $$
@@ -636,24 +640,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

@@ -21,6 +21,7 @@
 from process.data_structure.numerics import FiguresOfMerit, PROCESSRunMode
 from process.data_structure.physics_variables import DivertorNumberModels
 from process.data_structure.scan_variables import init_scan_variables
+from process.models.physics.profiles import DensityProfilePedestalType
 from process.models.stellarator.initialization import st_init
 from process.models.superconductors import (
     SuperconductorMaterial,
@@ -421,41 +422,49 @@ 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.physics.i_nd_plasma_pedestal_separatrix
+        )
         if (
-            data.physics.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.physics.nd_plasma_pedestal_electron
+            < data.physics.nd_plasma_separatrix_electron
+        ) or (
+            pedestal_type == DensityProfilePedestalType.GREENWALD_FRACTION
+            and data.physics.f_nd_plasma_pedestal_greenwald
+            < data.physics.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.physics.nd_plasma_pedestal_electron
-                < data.physics.nd_plasma_separatrix_electron
-            ):
-                raise ProcessValidationError(
-                    "Density pedestal is lower than separatrix density",
-                    nd_plasma_pedestal_electron=data.physics.nd_plasma_pedestal_electron,
-                    nd_plasma_separatrix_electron=data.physics.nd_plasma_separatrix_electron,
-                )
+            raise ProcessValidationError(
+                "Density pedestal is lower than separatrix density",
+                **(
+                    {
+                        "nd_plasma_pedestal_electron": data.physics.nd_plasma_pedestal_electron,
+                        "nd_plasma_separatrix_electron": data.physics.nd_plasma_separatrix_electron,
+                    }
+                    if pedestal_type == DensityProfilePedestalType.USER_INPUT
+                    else {
+                        "f_nd_plasma_pedestal_greenwald": data.physics.f_nd_plasma_pedestal_greenwald,
+                        "f_nd_plasma_separatrix_greenwald": data.physics.f_nd_plasma_separatrix_greenwald,
+                    }
+                ),
+            )
 
-            # Issue #589 Pedestal density is set manually using nd_plasma_pedestal_electron,
-            # but pedestal width = 0.
-            if (
-                abs(data.physics.radius_plasma_pedestal_density_norm - 1.0) <= 1e-7
-                and (
-                    data.physics.nd_plasma_pedestal_electron
-                    - data.physics.nd_plasma_separatrix_electron
-                )
-                >= 1e-7
-            ):
-                warn(
-                    "Density pedestal is at plasma edge "
-                    f"({data.physics.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron "
-                    f"({data.physics.nd_plasma_pedestal_electron}) differs from "
-                    f"nd_plasma_separatrix_electron ({data.physics.nd_plasma_separatrix_electron})",
-                    stacklevel=2,
-                )
+        if (
+            abs(data.physics.radius_plasma_pedestal_density_norm - 1.0) <= 1e-7
+            and (
+                data.physics.nd_plasma_pedestal_electron
+                - data.physics.nd_plasma_separatrix_electron
+            )
+            >= 1e-7
+        ):
+            warn(
+                "Density pedestal is at plasma edge "
+                f"({data.physics.radius_plasma_pedestal_density_norm = }), but nd_plasma_pedestal_electron "
+                f"({data.physics.nd_plasma_pedestal_electron}) differs from "
+                f"nd_plasma_separatrix_electron ({data.physics.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

@@ -149,9 +149,9 @@ def __post_init__(self):
     "f_p_div_lower": InputVariable("physics", float, range=(0.0, 1.0)),
     "f_plasma_fuel_deuterium": InputVariable("physics", float, range=(0.0, 1.0)),
     "ffwal": InputVariable("physics", float, range=(0.0, 10.0)),
-    "f_nd_plasma_pedestal_greenwald": InputVariable("physics", float, range=(-1.0, 5.0)),
+    "f_nd_plasma_pedestal_greenwald": InputVariable("physics", float, range=(0.1, 1.5)),
     "f_nd_plasma_separatrix_greenwald": InputVariable(
-        "physics", float, range=(-1.0, 5.0)
+        "physics", float, range=(0.001, 0.9)
     ),
     "f_plasma_fuel_helium3": InputVariable("physics", float, range=(-1.0, 5.0)),
     # TODO: does f_nd_impurity_electrons require an additional range?
@@ -578,6 +578,9 @@ def __post_init__(self):
     ),
     "nd_plasma_pedestal_electron": InputVariable("physics", float, range=(0.0, 1e21)),
     "nd_plasma_separatrix_electron": InputVariable("physics", 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)),
     "j_tf_coil_full_area": InputVariable("tfcoil", float, range=(10000.0, 1000000000.0)),
     "f_a_cs_turn_steel": InputVariable("pf_coil", float, range=(0.001, 0.999)),

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

@@ -3172,7 +3172,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(
@@ -3866,8 +3866,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([
@@ -4057,7 +4061,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}"
@@ -4068,7 +4075,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}"
         ),
@@ -4078,12 +4085,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,
@@ -4302,33 +4309,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

@@ -210,8 +210,8 @@ class IterationVariable:
         1.0e17,
         1.0e20,
     ),
-    145: IterationVariable("f_nd_plasma_pedestal_greenwald", "physics", 0.1, 0.9),
-    152: IterationVariable("f_nd_plasma_separatrix_greenwald", "physics", 0.001, 0.5),
+    145: IterationVariable("f_nd_plasma_pedestal_greenwald", "physics", 0.1, 1.5),
+    152: IterationVariable("f_nd_plasma_separatrix_greenwald", "physics", 0.001, 0.9),
     155: IterationVariable("pfusife", "ife", 5.0e2, 3.0e3),
     156: IterationVariable("rrin", "ife", 1.0, 1.0e1),
     158: IterationVariable(

process/data_structure/physics_variables.py

@@ -336,16 +336,16 @@ class PhysicsData:
     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 = 0.85
-    """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 = 0.50
-    """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`)
+    f_nd_plasma_separatrix_greenwald: float = 0.5
+    """Greenwald fraction of the separatrix density
     """
 
     f_plasma_fuel_helium3: float = 0.0
@@ -488,16 +488,22 @@ class PhysicsData:
     """
 
     nd_plasma_pedestal_electron: float = 4.0e19
-    """electron density of pedestal [m-3] (`i_plasma_pedestal==1)"""
+    """electron density of pedestal [m⁻³] (`i_plasma_pedestal==1)"""
 
     nd_plasma_separatrix_electron: float = 3.0e19
-    """electron density at separatrix [m-3] (`i_plasma_pedestal==1)"""
+    """electron density at separatrix [m⁻³] (`i_plasma_pedestal==1)"""
+
+    i_nd_plasma_pedestal_separatrix: int = 1
+    """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 = 0.0
     """critical ballooning parameter value"""
 
     nd_plasma_separatrix_electron_eich_max: float = 0.0
-    """Eich critical electron density at separatrix [m-3]"""
+    """Eich critical electron density at separatrix [m⁻³]"""
 
     plasma_res_factor: float = 1.0
     """plasma resistivity pre-factor"""
@@ -983,6 +989,9 @@ class PhysicsData:
     temp_plasma_electron_density_weighted_kev: float = 0.0
     """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 = 12.9
     """volume averaged ion temperature (keV). N.B. calculated from temp_plasma_electron_vol_avg_kev if `f_temp_plasma_ion_electron > 0.0`"""