Add enum for profile shape

process/models/physics/physics.py

@@ -35,6 +35,7 @@
 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
 
 if TYPE_CHECKING:
     from process.models.physics.plasma_current import (
@@ -336,19 +337,21 @@ def run(self):
         # If physics_variables.i_plasma_pedestal = 1 then set pedestal density to
         #   physics_variables.f_nd_plasma_pedestal_greenwald * Greenwald density limit
         # Note: this used to be done before plasma current
-        if (physics_variables.i_plasma_pedestal == 1) and (
-            physics_variables.f_nd_plasma_pedestal_greenwald >= 0e0
-        ):
+        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 (physics_variables.i_plasma_pedestal == 1) and (
-            physics_variables.f_nd_plasma_separatrix_greenwald >= 0e0
-        ):
+        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
@@ -2392,7 +2395,10 @@ def outplas(self):
             physics_variables.i_plasma_pedestal,
         )
 
-        if physics_variables.i_plasma_pedestal >= 1:
+        if (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PEDESTAL_PROFILE
+        ):
             if (
                 physics_variables.nd_plasma_electron_on_axis
                 < physics_variables.nd_plasma_pedestal_electron

process/models/physics/plasma_profiles.py

@@ -7,6 +7,7 @@
 from process.core import constants
 from process.core.exceptions import ProcessValueError
 from process.data_structure import divertor_variables, physics_variables
+from process.models.physics.profiles import PlasmaProfileShapeType
 
 logger = logging.getLogger(__name__)
 
@@ -59,7 +60,10 @@ def parameterise_plasma(self):
             )
 
         # Parabolic profile case
-        if physics_variables.i_plasma_pedestal == 0:
+        if (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PARABOLIC_PROFILE
+        ):
             self.parabolic_paramterisation()
             self.calculate_profile_factors()
             self.calculate_parabolic_profile_factors()
@@ -297,7 +301,10 @@ def calculate_parabolic_profile_factors():
         The function uses analytical parametric formulas to calculate the gradient information.
         The maximum normalized radius (rho_max) is obtained by equating the second derivative to zero.
         """
-        if physics_variables.i_plasma_pedestal == 0:
+        if (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PARABOLIC_PROFILE
+        ):
             if physics_variables.alphat > 1.0:
                 # Rho (normalized radius), where temperature derivative is largest
                 rho_te_max = 1.0 / np.sqrt(-1.0 + 2.0 * physics_variables.alphat)

process/models/physics/profiles.py

@@ -1,5 +1,6 @@
 import logging
 from abc import ABC, abstractmethod
+from enum import IntEnum
 
 import numpy as np
 import scipy as sp
@@ -9,6 +10,13 @@
 logger = logging.getLogger(__name__)
 
 
+class PlasmaProfileShapeType(IntEnum):
+    """Enum for i_plasma_pedestal method types"""
+
+    PARABOLIC_PROFILE = 0
+    PEDESTAL_PROFILE = 1
+
+
 class Profile(ABC):
     """Abstract base class used to create and hold profiles (temperature, density)"""
 
@@ -123,7 +131,10 @@ def calculate_profile_y(
             Density peaking parameter.
         """
 
-        if physics_variables.i_plasma_pedestal == 0:
+        if (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PARABOLIC_PROFILE
+        ):
             self.profile_y = n0 * (1 - rho**2) ** alphan
 
         # Input checks
@@ -212,12 +223,18 @@ def ncore(
     def set_physics_variables(self):
         """Calculates and sets physics variables required for the profile."""
 
-        if physics_variables.i_plasma_pedestal == 0:
+        if (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PARABOLIC_PROFILE
+        ):
             physics_variables.nd_plasma_electron_on_axis = (
                 physics_variables.nd_plasma_electrons_vol_avg
                 * (1.0 + physics_variables.alphan)
             )
-        elif physics_variables.i_plasma_pedestal == 1:
+        elif (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PEDESTAL_PROFILE
+        ):
             physics_variables.nd_plasma_electron_on_axis = self.ncore(
                 physics_variables.radius_plasma_pedestal_density_norm,
                 physics_variables.nd_plasma_pedestal_electron,
@@ -284,7 +301,10 @@ def calculate_profile_y(
         References:
             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
         """
-        if physics_variables.i_plasma_pedestal == 0:
+        if (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PARABOLIC_PROFILE
+        ):
             # profile values of 0 cause divide by 0 errors so ensure the profile value is at least 1e-8
             # which is small enough that it won't make a difference to any calculations
             self.profile_y = np.maximum(t0 * (1 - rho**2) ** alphat, 1e-8)
@@ -374,12 +394,18 @@ def tcore(
 
     def set_physics_variables(self):
         """Calculates and sets physics variables required for the temperature profile."""
-        if physics_variables.i_plasma_pedestal == 0:
+        if (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PARABOLIC_PROFILE
+        ):
             physics_variables.temp_plasma_electron_on_axis_kev = (
                 physics_variables.temp_plasma_electron_vol_avg_kev
                 * (1.0 + physics_variables.alphat)
             )
-        elif physics_variables.i_plasma_pedestal == 1:
+        elif (
+            PlasmaProfileShapeType(physics_variables.i_plasma_pedestal)
+            == PlasmaProfileShapeType.PEDESTAL_PROFILE
+        ):
             physics_variables.temp_plasma_electron_on_axis_kev = self.tcore(
                 physics_variables.radius_plasma_pedestal_temp_norm,
                 physics_variables.temp_plasma_pedestal_kev,