🚧 3407 update the density limit section of the docs to be up to date and include all models

documentation/proc-pages/physics-models/plasma_density.md

@@ -1,19 +1,125 @@
 # Density Limit
 
-Several density limit models[^1] are available in PROCESS. These are
-calculated in routine `culdlm`, which is called by `physics`. To enforce any of 
-these limits, turn on constraint equation no. 5 with iteration variable no. 9 
-(`fdene`). In addition, switch `idensl` must be set to the relevant value, as 
-follows:
-
-| `idensl` | Description |
-| :-: | - |
-| 1 | ASDEX model |
-| 2 | Borrass model for ITER, I |
-| 3 | Borrass model for ITER, II |
-| 4 | JET edge radiation model |
-| 5 | JET simplified model |
-| 6 | Hugill-Murakami $M.q$ model |
-| 7 | Greenwald model: $n_G=10^{14} \frac{I_p}{\pi a^2}$ where the units are m and ampere. For the Greenwald model the limit applies to the line-averaged electron density, not the volume-averaged density. |
-
-[^1]: T. C. Hender et al., 'Physics Assessment for the European Reactor Study',
+Several density limit models are available in PROCESS. These are
+calculated in routine `calculate_density_limit()`, which is called by `physics`.
+
+This constraint can be activated by stating `icc = 5` in the input file.
+
+The value of `i_density_limit` can be set to apply the relevant limit . The scaling value `fdene` can be varied also. 
+
+For the `i_density_limit = 1-5,8` scalings we scale the function output by the separatrix to volume averaged electron density so that we can set the limit on the volume averaged. **Therefore it is recommended to only use these scalings with an H-mode profile (`ipedestal == 1`) otherwise the separatrix density (`nesep`) will not be calculated.**
+
+For the models below $P_{\perp}$ is the mean heat flux density across the separatrix ($\mathrm{MW}/\mathrm{m^2}$), which we take as the divertor power divided by the plasma surface area.
+
+-----------------
+
+## ASDEX model
+
+Switch value: `i_density_limit = 1`[^1][^2]
+
+$$
+n_{\text{b}}^{\text{crit}} = 1.54 \frac{P_{\perp}^{0.43}B_{\text{T}}^{0.31}}{\left(q_{95}R\right)^{0.45}}
+$$
+
+-----------------
+
+## Borrass model for ITER, I
+
+Switch value: `i_density_limit = 2` [^1]
+
+$$
+n_{\text{b}}^{\text{crit}} = C \frac{P_{\perp}^{0.53}B_{\text{T}}^{0.31}}{\left(q_{95}R\right)^{0.22}}
+$$
+
+$C \approx$  1.8 for ITER-like conditions.
+
+-----------------
+
+## Borrass model for ITER, II 
+
+Switch value: `i_density_limit = 3` [^1]
+
+$$
+n_{\text{b}}^{\text{crit}} = 0.5 \frac{P_{\perp}^{0.57}B_{\text{T}}^{0.31}}{\left(q_{95}R\right)^{0.09}}
+$$
+
+-----------------
+
+## JET edge radiation model
+
+Switch value: `i_density_limit = 4` [^1]
+
+$$
+n_{\text{b}}^{\text{crit}} = P_{\text{in}}^{0.5} \frac{1}{\left[\left(Z_{\text{eff}}-1\right)\left(1-\frac{4}{3q_{\text{c}}}\right)\right]^{0.5}}
+$$
+
+-----------------
+
+## JET simplified model
+
+Switch value: `i_density_limit = 5` [^1]
+
+$$
+n_{\text{b}}^{\text{crit}} = 0.237 P^{0.5}
+$$
+
+For a radiation from a shell thickness $\Delta$, this may be written as:
+
+$$
+n_{\text{b}}^{\text{crit}} = 0.147 \frac{P^{0.5}}{\left[Ra\Delta \sqrt{\frac{\left(1+\kappa \right)}{2}}\right]^{0.5}}
+$$
+
+where $\kappa \approx 1.5, \Delta \approx 0.1a$ has been taken from JET.
+
+-----------------
+
+## Hugill-Murakami model
+
+Switch value: `i_density_limit = 6` [^2]
+
+$$
+\langle n_{^{\text{crit}}} \rangle \approx \frac{3.0 B_{\text{T}}}{R_0 q_{\text{cyl}}}
+$$
+
+
+-----------------
+
+## Greenwald model
+
+Switch value: `i_density_limit = 7` [^3][^4]
+
+$$
+\overline{n}_{\text{e}}^{\text{ crit}} = 1.0 \times 10^{14} \frac{I_\text{p}}{\pi a^2}
+$$
+
+For the Greenwald model the limit applies to the line-averaged electron density, not the volume-averaged density. The plasma current term is given in $[\mathrm{A}]$ and the minor radius in $[\mathrm{m}]$
+
+---------------------
+
+## ASDEX New model
+
+Switch value: `i_density_limit = 8` [^5][^6]
+
+$$
+\overline{n}_{\text{sep}}^{\text{ crit}} = 1.0 \times 10^{20} \times 0.506 \pm 0.192 \frac{P_\text{heat}^{0.396\pm0.13} I_{\text{p}}^{0.265\pm 0.14}}{q_{95}^{0.323 \pm 0.14}}
+$$
+
+-----------------
+
+## Key Constraints
+
+[^1]: T.C.Hender et.al., 'Physics Assesment of the European Reactor Study', AEA FUS 172, 1992
+
+[^2]: N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989',
+
+[^3]: M. Greenwald et al., “A new look at density limits in tokamaks,” Nuclear Fusion, vol. 28, no. 12, pp. 2199–2207, Dec. 1988, doi: https://doi.org/10.1088/0029-5515/28/12/009.
+
+[^4]: M. Greenwald, “Density limits in toroidal plasmas,” Plasma Physics and Controlled Fusion, vol. 44, no. 8, pp. R27–R53, Jul. 2002, doi: https://doi.org/10.1088/0741-3335/44/8/201.
+
+[^5]: J. W. Berkery et al., “Density limits as disruption forecasters for spherical tokamaks,” Plasma Physics and Controlled Fusion, vol. 65, no. 9, pp. 095003–095003, Jul. 2023, doi: https://doi.org/10.1088/1361-6587/ace476.
+
+[^6]: M. Bernert et al., “The H-mode density limit in the full tungsten ASDEX Upgrade tokamak,” vol. 57, no. 1, pp. 014038–014038, Nov. 2014, doi: https://doi.org/10.1088/0741-3335/57/1/014038.
+‌
+‌
+‌
+‌

examples/data/csv_output_large_tokamak_MFILE.DAT

@@ -1577,7 +1577,7 @@ iculbl = 1
 i_plasma_current    = 4
 
 * Switch for density limit to enforce
-idensl = 7
+i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
 ifalphap = 1

examples/data/large_tokamak_1_MFILE.DAT

@@ -1571,7 +1571,7 @@ iculbl = 1
 i_plasma_current    = 4
 
 * Switch for density limit to enforce
-idensl = 7
+i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
 ifalphap = 1

examples/data/large_tokamak_2_MFILE.DAT

@@ -1571,7 +1571,7 @@ iculbl = 1
 i_plasma_current    = 4
 
 * Switch for density limit to enforce
-idensl = 7
+i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
 ifalphap = 1

examples/data/large_tokamak_3_MFILE.DAT

@@ -1571,7 +1571,7 @@ iculbl = 1
 i_plasma_current    = 4
 
 * Switch for density limit to enforce
-idensl = 7
+i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
 ifalphap = 1

examples/data/large_tokamak_4_MFILE.DAT

@@ -1571,7 +1571,7 @@ iculbl = 1
 i_plasma_current    = 4
 
 * Switch for density limit to enforce
-idensl = 7
+i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
 ifalphap = 1

examples/data/large_tokamak_IN.DAT

@@ -381,7 +381,7 @@ iculbl = 1
 i_plasma_current    = 4
 
 * Switch for density limit to enforce
-idensl = 7
+i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
 ifalphap = 1

examples/data/scan_MFILE.DAT

@@ -9241,7 +9241,7 @@ hfact    = 1.1 * H factor on energy confinement times (iteration variable 10)
 i_bootstrap_current     = 4 * Switch for bootstrap current scaling;
 iculbl   = 1 * Switch for beta limit scaling (constraint equation 24);
 i_plasma_current    = 4 * Switch for plasma current scaling to use;
-idensl   = 7 * Switch for density limit to enforce (constraint equation 5);
+i_density_limit   = 7 * Switch for density limit to enforce (constraint equation 5);
 ifalphap = 1 * Switch for fast alpha pressure calculation;
 ifispact = 0 * Switch for neutronics calculations;
 iinvqd   = 1 * Switch for inverse quadrature in l-mode scaling laws 5 and 9;

examples/data/scan_example_file_IN.DAT

@@ -381,7 +381,7 @@ iculbl = 1
 i_plasma_current    = 4
 
 * Switch for density limit to enforce
-idensl = 7
+i_density_limit = 7
 
 * Switch for fast alpha pressure calculation
 ifalphap = 1

process/io/obsolete_vars.py

@@ -127,6 +127,7 @@
     "ftrit": "f_tritium",
     "fhe3": "f_helium3",
     "falpha": "f_alpha_plasma",
+    "idensl": "i_density_limit",
     "ftburn": "ft_burn",
     "ftohs": "ft_current_ramp_up",
     "tbrnmn": "t_burn_min",

process/io/plot_proc.py

@@ -3110,11 +3110,96 @@ def plot_h_threshold_comparison(
     axis.set_facecolor("#f0f0f0")
 
 
+def plot_density_limit_comparison(
+    axis: plt.Axes, mfile_data: mf.MFile, scan: int
+) -> None:
+    """
+    Function to plot a scatter box plot of different density limit comparisons.
+
+    Arguments:
+        axis (plt.Axes): Axis object to plot to.
+        mfile_data (mf.MFile): MFILE data object.
+        scan (int): Scan number to use.
+    """
+    old_asdex = mfile_data.data["dlimit(1)"].get_scan(scan)
+    borrass_iter_i = mfile_data.data["dlimit(2)"].get_scan(scan)
+    borrass_iter_ii = mfile_data.data["dlimit(3)"].get_scan(scan)
+    jet_edge_radiation = mfile_data.data["dlimit(4)"].get_scan(scan)
+    jet_simplified = mfile_data.data["dlimit(5)"].get_scan(scan)
+    hugill_murakami = mfile_data.data["dlimit(6)"].get_scan(scan)
+    greenwald = mfile_data.data["dlimit(7)"].get_scan(scan)
+    asdex_new = mfile_data.data["dlimit(8)"].get_scan(scan)
+
+    # Data for the box plot
+    data = {
+        "Old ASDEX": old_asdex,
+        "Borrass ITER I": borrass_iter_i,
+        "Borrass ITER II": borrass_iter_ii,
+        "JET Edge Radiation": jet_edge_radiation,
+        "JET Simplified": jet_simplified,
+        "Hugill-Murakami": hugill_murakami,
+        "Greenwald": greenwald,
+        "ASDEX New": asdex_new,
+    }
+
+    # 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_density_limit = np.mean(data_values)
+    std_density_limit = np.std(data_values)
+    median_density_limit = np.median(data_values)
+
+    # Plot average, standard deviation, and median as text
+    axis.text(
+        1.02,
+        0.2,
+        rf"Average: {avg_density_limit * 1e-20:.4f} $\times 10^{{20}}$",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+    axis.text(
+        1.02,
+        0.15,
+        rf"Standard Dev: {std_density_limit * 1e-20:.4f} $\times 10^{{20}}$",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+    axis.text(
+        1.02,
+        0.1,
+        rf"Median: {median_density_limit * 1e-20:.4f} $\times 10^{{20}}$",
+        transform=axis.transAxes,
+        fontsize=9,
+    )
+
+    axis.set_title("Density Limit Comparison")
+    axis.set_ylabel(r"Density Limit [$10^{20}$ m$^{-3}$]")
+    axis.yaxis.set_major_formatter(plt.FuncFormatter(lambda x, _: f"{x * 1e-20:.1f}"))
+    axis.set_xlim([0.5, 1.5])
+    axis.set_xticks([])
+    axis.set_xticklabels([])
+    axis.set_facecolor("#f0f0f0")
+
+
 def main_plot(
     fig1,
     fig2,
     fig3,
     fig4,
+    fig5,
     m_file_data,
     scan,
     imp="../data/lz_non_corona_14_elements/",
@@ -3217,6 +3302,9 @@ def main_plot(
     plot_10 = fig4.add_subplot(224)
     plot_h_threshold_comparison(plot_10, m_file_data, scan)
 
+    plot_11 = fig5.add_subplot(221)
+    plot_density_limit_comparison(plot_11, m_file_data, scan)
+
 
 def main(args=None):
     # TODO The use of globals here isn't ideal, but is required to get main()
@@ -3472,13 +3560,15 @@ def main(args=None):
     page2 = plt.figure(figsize=(12, 9), dpi=80)
     page3 = plt.figure(figsize=(12, 9), dpi=80)
     page4 = plt.figure(figsize=(12, 9), dpi=80)
+    page5 = plt.figure(figsize=(12, 9), dpi=80)
 
     # run main_plot
     main_plot(
         page1,
         page2,
         page3,
         page4,
+        page5,
         m_file,
         scan=scan,
         demo_ranges=demo_ranges,
@@ -3491,6 +3581,7 @@ def main(args=None):
         pdf.savefig(page2)
         pdf.savefig(page3)
         pdf.savefig(page4)
+        pdf.savefig(page5)
 
     # show fig if option used
     if args.show:
@@ -3500,6 +3591,7 @@ def main(args=None):
     plt.close(page2)
     plt.close(page3)
     plt.close(page4)
+    plt.close(page5)
 
 
 if __name__ == "__main__":