Add SOL physics class

documentation/source/physics-models/plasma_scrape_off_layer.md

@@ -0,0 +1,60 @@
+# Scrape off Layer | `ScrapeOffLayer`
+
+## Power Decay Lengths
+
+In tokamak scrape-off layer (SOL) physics, $\lambda_q$ is the radial heat flux e-folding length, representing the narrow width over which exhaust power drops exponentially. It measures how tightly thermal energy exiting the core plasma is compressed onto open magnetic field lines.
+
+-----------
+
+### Eich 2013 Model | `calculate_eich2013_sol_power_decay_length()`
+
+The power decay length in metres is given by[^eich_2013]:
+
+$$
+\lambda_q = 1.35 \times 10^{-3} P_{\text{sep}}^{-0.02}R_0^{0.04}B_{\text{p}}(a)^{-0.92}\epsilon^{0.42}
+$$
+
+Here $P_{\text{sep}}$ is the plasma separatrix power in $\text{MW}$, $R_0$ is the plasma major radius in $\text{m}$, $B_{\text{p}}(a)$ is the plasma poloidal field at the outboard mid-plane measured in $\text{T}$ and $\epsilon$ is the plasma inverse aspect ratio.
+
+This can be found in Table 3 from Eich et.al [^eich_2013]
+
+The $R^2$ value for this fit is 0.88
+
+----------
+
+### MAST 2014 I | `calculate_mast2014_sol_power_decay_length_1()`
+
+The power decay length in metres is given by[^mast_2014]:
+
+$$
+\lambda_q = 1.84(\pm0.48) \times 10^{-3} P_{\text{sep}}^{0.18(\pm0.07)}B_{\text{p}}(a)^{-0.68(\pm0.14)}
+$$
+
+Here $P_{\text{sep}}$ is the plasma separatrix power in $\text{MW}$, and $B_{\text{p}}(a)$ is the plasma poloidal field at the outboard mid-plane measured in $\text{T}$.
+
+This can be found in Table 2 and Equation 3 from Thornton et.al [^mast_2014]
+
+The $R^2$ value for this fit is 0.56
+
+----------
+
+### MAST 2014 II | `calculate_mast2014_sol_power_decay_length_2()`
+
+The power decay length in metres is given by[^mast_2014]:
+
+$$
+\lambda_q = 4.57(\pm0.54) \times 10^{-3} P_{\text{sep}}^{0.22(\pm0.08)}I_{\text{p}}^{-0.64(\pm0.15)}
+$$
+
+Here $P_{\text{sep}}$ is the plasma separatrix power in $\text{MW}$, and $I_{\text{p}}$ is the total plasma current measured in $\text{MA}$.
+
+This can be found in Table 2 and Equation 4 from Thornton et.al [^mast_2014]
+
+The $R^2$ value for this fit is 0.55
+
+--------
+
+[^eich_2013]: T. Eich et al., “Scaling of the tokamak near the scrape-off layer H-mode power width and implications for ITER,” Nuclear Fusion, vol. 53, no. 9 p. 093031, Aug. 2013, doi: 10.1088/0029-5515/53/9/093031.
+
+[^mast_2014]: A. J. Thornton and A. Kirk, “Scaling of the scrape-off layer width during inter-ELM H modes on MAST as measured by infrared thermography,”
+Plasma Physics and Controlled Fusion, vol. 56, no. 5, p. 055008, Apr. 2014, doi: 10.1088/0741-3335/56/5/055008.

mkdocs.yml

@@ -67,6 +67,7 @@ nav:
               - Confinement time: physics-models/plasma_confinement.md
               - L-H transition: physics-models/plasma_h_mode.md
               - Plasma Core Power Balance: physics-models/plasma_power_balance.md
+              - Plasma Scrape-off Layer: physics-models/plasma_scrape_off_layer.md
               - Detailed Plasma Physics: physics-models/detailed_physics.md
           - Pulsed Plant Operation: physics-models/pulsed-plant.md
       - Engineering Models:

process/core/io/plot/summary.py

@@ -8914,6 +8914,85 @@ def plot_bootstrap_comparison(axis: plt.Axes, mfile: MFile, scan: int):
     axis.set_facecolor("#f0f0f0")
 
 
+def plot_sol_power_decay_length_comparison(axis: plt.Axes, mfile: MFile, scan: int):
+    """Function to plot a scatter box plot of SOL power decay lengths (λ_q).
+
+    Parameters
+    ----------
+    axis :
+        axis object to plot to
+    mfile :
+        MFILE data object
+    scan :
+        scan number to use
+    """
+    len_plasma_sol_eich13_power_decay_mm = (
+        mfile.get("len_plasma_sol_eich13_power_decay", scan=scan) * 1e3
+    )
+    len_plasma_sol_mast14_power_decay_1_mm = (
+        mfile.get("len_plasma_sol_mast14_power_decay_1", scan=scan) * 1e3
+    )
+    len_plasma_sol_mast14_power_decay_2_mm = (
+        mfile.get("len_plasma_sol_mast14_power_decay_2", scan=scan) * 1e3
+    )
+
+    # Data for the box plot
+    data = {
+        "Eich 2013": len_plasma_sol_eich13_power_decay_mm,
+        "MAST 2014 (1)": len_plasma_sol_mast14_power_decay_1_mm,
+        "MAST 2014 (2)": len_plasma_sol_mast14_power_decay_2_mm,
+    }
+    # Create the violin plot
+    axis.violinplot(data.values(), showextrema=False)
+
+    # Create the box plot
+    axis.boxplot(
+        data.values(), showfliers=True, showmeans=True, meanline=True, widths=0.3
+    )
+
+    # Scatter plot for each data point
+    colors = plt.cm.plasma(np.linspace(0, 1, len(data.values())))
+    for index, (key, value) in enumerate(data.items()):
+        axis.scatter(1, value, color=colors[index], label=key, alpha=1.0)
+    axis.legend(loc="upper left", bbox_to_anchor=(1, 1))
+
+    # Calculate average, standard deviation, and median
+    data_values = list(data.values())
+    avg_decay_length = np.mean(data_values)
+    std_decay_length = np.std(data_values)
+    median_decay_length = np.median(data_values)
+
+    # Plot average, standard deviation, and median as text
+    axis.text(
+        1.02,
+        0.2,
+        f"Average: {avg_decay_length:.4f}",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+    axis.text(
+        1.02,
+        0.15,
+        f"Standard Dev: {std_decay_length:.4f}",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+    axis.text(
+        1.02,
+        0.1,
+        f"Median: {median_decay_length:.4f}",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+
+    axis.set_title("SOL Power Decay Length ($\\lambda_q$) Comparison")
+    axis.set_ylabel("Power Decay Length [mm]")
+    axis.set_xlim([0.5, 1.5])
+    axis.set_xticks([])
+    axis.set_xticklabels([])
+    axis.set_facecolor("#f0f0f0")
+
+
 def plot_h_threshold_comparison(axis: plt.Axes, mfile: MFile, scan: int, u_seed=None):
     """Function to plot a scatter box plot of L-H threshold power comparisons.
 
@@ -16090,6 +16169,10 @@ def _add_page(name: str | None = None):
         pages["plasma_compare_2"].add_subplot(224), m_file, scan
     )
 
+    plot_sol_power_decay_length_comparison(
+        _add_page("plasma_compare_3").add_subplot(221), m_file, scan
+    )
+
     plot_debye_length_profile(
         _add_page("microscopic_quantities").add_subplot(232), m_file, scan
     )

process/data_structure/physics_variables.py

@@ -1651,6 +1651,15 @@ class PhysicsData:
     res_plasma_fuel_spitzer_vol_avg: float = 0.0
     """Volume averaged plasma Spitzer resistivity due to fuel ions (ohm m)"""
 
+    len_plasma_sol_eich13_power_decay: float = 0.0
+    """Eich 2013 power decay length in the scrape-off layer scaling (λ_q) [m]"""
+
+    len_plasma_sol_mast14_power_decay_1: float = 0.0
+    """MAST 2014 power decay length in the scrape-off layer scaling 1 (λ_q) [m]"""
+
+    len_plasma_sol_mast14_power_decay_2: float = 0.0
+    """MAST 2014 power decay length in the scrape-off layer scaling 2 (λ_q) [m]"""
+
     dt_power_density_plasma: float = 0.0
     sigmav_dt_average: float = 0.0
     dhe3_power_density: float = 0.0

process/main.py

@@ -103,6 +103,7 @@
 from process.models.physics.plasma_geometry import PlasmaGeom
 from process.models.physics.plasma_profiles import PlasmaProfile
 from process.models.physics.profiles import NeProfile, TeProfile
+from process.models.physics.scrape_off_layer import ScrapeOffLayer
 from process.models.power import Power
 from process.models.pulse import Pulse
 from process.models.shield import Shield
@@ -679,6 +680,7 @@ def __init__(self, data: DataStructure):
         self.plasma_current = PlasmaCurrent()
         self.plasma_fields = PlasmaFields()
         self.plasma_dia_current = PlasmaDiamagneticCurrent()
+        self.scrape_off_layer = ScrapeOffLayer()
         self.physics = Physics(
             plasma_profile=self.plasma_profile,
             current_drive=self.current_drive,
@@ -693,6 +695,7 @@ def __init__(self, data: DataStructure):
             plasma_fields=self.plasma_fields,
             plasma_dia_current=self.plasma_dia_current,
             plasma_geometry=self.plasma_geom,
+            scrape_off_layer=self.scrape_off_layer,
         )
         self.physics_detailed = DetailedPhysics(
             plasma_profile=self.plasma_profile,
@@ -782,6 +785,7 @@ def models(self) -> tuple[Model, ...]:
             self.plasma_bootstrap_current,
             self.plasma_exhaust,
             self.plasma_current,
+            self.scrape_off_layer,
             self.neoclassics,
             self.plasma_inductance,
             self.ne_profile,

process/models/physics/physics.py

@@ -45,6 +45,7 @@
     from process.models.physics.plasma_fields import PlasmaFields
     from process.models.physics.plasma_geometry import PlasmaGeom
     from process.models.physics.plasma_profiles import PlasmaProfile
+    from process.models.physics.scrape_off_layer import ScrapeOffLayer
 
 logger = logging.getLogger(__name__)
 
@@ -195,6 +196,7 @@ def __init__(
         plasma_fields: PlasmaFields,
         plasma_dia_current: PlasmaDiamagneticCurrent,
         plasma_geometry: PlasmaGeom,
+        scrape_off_layer: ScrapeOffLayer,
     ):
         self.outfile = constants.NOUT
         self.mfile = constants.MFILE
@@ -211,6 +213,7 @@ def __init__(
         self.fields = plasma_fields
         self.dia_current = plasma_dia_current
         self.geometry = plasma_geometry
+        self.scrape_off_layer = scrape_off_layer
 
     def output(self) -> None:
         """Output plasma physics information."""
@@ -813,6 +816,13 @@ def run(self):
             )
         )
 
+        # ===============================
+
+        # Calculate scrape-off layer physics
+
+        self.scrape_off_layer.run()
+
+        # ===============================
         self.data.physics.pflux_plasma_surface_neutron_avg_mw = (
             self.data.physics.p_neutron_total_mw / self.data.physics.a_plasma_surface
         )
@@ -2334,6 +2344,8 @@ def outplas(self):
 
         self.exhaust.output()
 
+        self.scrape_off_layer.output()
+
         if self.data.stellarator.istell == 0:
             self.plasma_transition.output()