3145 sauter geometry without using the sauter plasma current

documentation/proc-pages/eng-models/divertor.md

@@ -5,9 +5,9 @@ The principal outputs from the code are the divertor heat load, used to
 determine its lifetime, and its peak temperature. The divertor is cooled either 
 by gaseous helium or by pressurised water.
 
-Switch `snull` controls the overall plasma configuration. Setting `snull = 0` 
+Switch `i_single_null` controls the overall plasma configuration. Setting `i_single_null= 0` 
 corresponds to an up-down symmetric, double null configuration, while 
-`snull = 1` (the default) assumes a single null plasma with the divertor at the 
+`i_single_null= 1` (the default) assumes a single null plasma with the divertor at the 
 bottom of the machine. The vertical build and PF coil current scaling 
 algorithms take the value of this switch into account, although not the plasma 
 geometry at present.

documentation/proc-pages/eng-models/pf-coil.md

@@ -25,7 +25,7 @@ values for `rpf1`, `rpf2`, `zref(j)` and `routr` should be adjusted by the user
 coils accurately.
 
 The three possible values of `ipfloc(j)` correspond to the following PF coil positions: (Redo taking 
-into account snull and other recent changes e.g. rclsnorm)
+into account `i_single_null` and other recent changes e.g. rclsnorm)
 
 `ipfloc(j)` = 1: PF coils are placed above the central solenoid (one group only);
 *R* = `rohc` + `rpf1`<br>

documentation/proc-pages/fusion-devices/spherical-tokamak.md

@@ -10,7 +10,7 @@
         title="Schematic diagram of the Power Core radial build" 
         width="650" height="100" />
         <br><br>
-        <figcaption><i>Figure 1: Schematic diagram of the fusion power core of a typical tokamak power plant modelled by `PROCESS`, showing the relative positions of the components. A double null plasma is assumed (`snull=0`) - compare Figure 2, and the first wall, blanket, shield and vacuum vessel are D-shaped in cross-section (chosen by setting switch `fwbsshape=1`) - compare Figure 3. Also shown are the code variables used to define the thicknesses of the components. The arrowed labels adjacent to the axes are the total 'builds' to that point. The precise locations and sizes of the PF coils are calculated within the code.
+        <figcaption><i>Figure 1: Schematic diagram of the fusion power core of a typical tokamak power plant modelled by `PROCESS`, showing the relative positions of the components. A double null plasma is assumed (`i_single_null=0`) - compare Figure 2, and the first wall, blanket, shield and vacuum vessel are D-shaped in cross-section (chosen by setting switch `fwbsshape=1`) - compare Figure 3. Also shown are the code variables used to define the thicknesses of the components. The arrowed labels adjacent to the axes are the total 'builds' to that point. The precise locations and sizes of the PF coils are calculated within the code.
         </i></figcaption>
         <br>
         </center>
@@ -42,7 +42,7 @@ Switch `itart` provides overall control of the ST switches within the code, and
 
 | Switch | Conventional aspect ratio (`itart = 0`) | Low aspect ratio (`itart = 1`) |
 | --- | --- | --- |
-| `ishape` | 0, 2, 3, 4 | 1, 5, 6, 7, 8 |
+| `i_plasma_geometry` | 0, 2, 3, 4 | 1, 5, 6, 7, 8 |
 | `i_bootstrap_current` | 1, 2, 3 | 2, 3 |
 | `i_plasma_current` | 1, 3, 4, 5, 6, 7 | 2, 9 |
 | `itfsup` | 0, 1 | 0 |
@@ -55,7 +55,7 @@ Switch `itart` provides overall control of the ST switches within the code, and
     title="Schematic diagram of the Power Core radial build" 
     width="650" height="100" />
     <br><br>
-    <figcaption><i>Figure 2: Schematic diagram of the fusion power core of a typical tokamak power plant modelled by `PROCESS`, showing the relative positions of the components. A single null plasma is assumed ('snull=1') - compare Figure 3. The radial build is the same as for a double null configuration; shown along the vertical axis are the code variables used to define the vertical thicknesses of the components. The arrow labels adjacent to the axis are the total 'builds' (distance from the midplane, Z=0) to that point. The precise locations and sizes of the PF coils are calculated within the code.
+    <figcaption><i>Figure 2: Schematic diagram of the fusion power core of a typical tokamak power plant modelled by `PROCESS`, showing the relative positions of the components. A single null plasma is assumed ('i_single_null=1') - compare Figure 3. The radial build is the same as for a double null configuration; shown along the vertical axis are the code variables used to define the vertical thicknesses of the components. The arrow labels adjacent to the axis are the total 'builds' (distance from the midplane, Z=0) to that point. The precise locations and sizes of the PF coils are calculated within the code.
     </i></figcaption>
     <br>
     </center>

documentation/proc-pages/physics-models/error.txt

@@ -119,15 +119,15 @@ values. The shape of the plasma cross-section is given by its last
 closed flux surface (LCFS) elongation \(\kappa\) (\texttt{kappa}) and
 triangularity \(\delta\) (\texttt{triang}), which can be scaled
 automatically with the aspect ratio if required using switch
-\texttt{ishape}:
+\texttt{i_plasma_geometry}:
 
 \begin{itemize}
 \tightlist
 \item
-  \texttt{ishape\ =\ 0} -- the input values for \texttt{kappa} and
+  \texttt{i_plasma_geometry\ =\ 0} -- the input values for \texttt{kappa} and
   \texttt{triang} are used directly.
 \item
-  \texttt{ishape\ =\ 1} -- the following scaling is used, which is
+  \texttt{i_plasma_geometry\ =\ 1} -- the following scaling is used, which is
   suitable for low aspect ratio machines (\(\epsilon = 1/A\)) \footnote{J.D.
     Galambos, `STAR Code : Spherical Tokamak Analysis and Reactor Code',
     Unpublished internal Oak Ridge document.}:
@@ -141,7 +141,7 @@ automatically with the aspect ratio if required using switch
 \begin{itemize}
 \tightlist
 \item
-  \texttt{ishape\ =\ 2} -- the Zohm ITER scaling \footnote{H. Zohm et
+  \texttt{i_plasma_geometry\ =\ 2} -- the Zohm ITER scaling \footnote{H. Zohm et
     al, `On the Physics Guidelines for a Tokamak DEMO', FTP/3-3, Proc.
     IAEA Fusion Energy Conference, October 2012, San Diego} is used to
   calculate the elongation: \[
@@ -151,7 +151,7 @@ automatically with the aspect ratio if required using switch
   unchanged.
 \end{itemize}
 
-If \texttt{ishape\ =\ 0,\ 1,\ 2}, the values for the plasma shaping
+If \texttt{i_plasma_geometry\ =\ 0,\ 1,\ 2}, the values for the plasma shaping
 parameters at the 95\% flux surface are calculated as follows:
 
 \[\begin{aligned}
@@ -162,13 +162,13 @@ parameters at the 95\% flux surface are calculated as follows:
 \begin{itemize}
 \tightlist
 \item
-  If \texttt{ishape\ =\ 3}, the Zohm ITER scaling is used to calculate
-  the elongation (as for \texttt{ishape\ =\ 2} above), but the
+  If \texttt{i_plasma_geometry\ =\ 3}, the Zohm ITER scaling is used to calculate
+  the elongation (as for \texttt{i_plasma_geometry\ =\ 2} above), but the
   triangularity at the 95\% flux surface is input via variable
   \texttt{triang95}, and the LCFS triangularity \texttt{triang} is
   calculated from it, rather than the other way round.
 \item
-  Finally, if \texttt{ishape\ =\ 4}, the 95\% flux surface values
+  Finally, if \texttt{i_plasma_geometry\ =\ 4}, the 95\% flux surface values
   \texttt{kappa95} and \texttt{triang95} are both used as inputs, and
   the LCFS values are calculated from them by inverting the equations
   above.
@@ -182,10 +182,10 @@ significant elongation \footnote{Y. Sakamoto, `Recent progress in
   vertical stability analysis in JA', Task meeting EU-JA \#16, Fusion
   for Energy, Garching, 24--25 June 2014}. The maximum distance
 \(r_{\text{shell, max}}\) of this shell from the centre of the plasma
-may be set using input parameter \texttt{cwrmax}, such that
-\(r_{\text{shell, max}} =\) \texttt{cwrmax*rminor}. Constraint equation
+may be set using input parameter \texttt{f_r_conducting_wall}, such that
+\(r_{\text{shell, max}} =\) \texttt{f_r_conducting_wall*rminor}. Constraint equation
 no. 23 should be turned on with iteration variable no.~104
-(\texttt{fcwr}) to enforce this. (A scaling of \texttt{cwrmax} with
+(\texttt{fcwr}) to enforce this. (A scaling of \texttt{f_r_conducting_wall} with
 elongation should be available shortly.)
 
 The plasma surface area, cross-sectional area and volume are calculated
@@ -994,18 +994,18 @@ is derived directly from the energy confinement scaling law.
 \texttt{iradloss\ =\ 0} -- Total power lost is scaling power plus
 radiation
 
-\texttt{pscaling\ +\ pradpv\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/plasma_volume}
+\texttt{pscaling\ +\ pradpv\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/vol_plasma}
 
 \texttt{iradloss\ =\ 1} -- Total power lost is scaling power plus core
 radiation only
 
-\texttt{pscaling\ +\ pcoreradpv\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/plasma_volume}
+\texttt{pscaling\ +\ pcoreradpv\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/vol_plasma}
 
 \texttt{iradloss\ =\ 2} -- Total power lost is scaling power only, with
 no additional allowance for radiation. This is not recommended for power
 plant models.
 
-\texttt{pscaling\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/plasma_volume}
+\texttt{pscaling\ =\ f_alpha_plasma*alpha_power_density\ +\ charged_power_density\ +\ pden_plasma_ohmic_mw\ +\ pinjmw/vol_plasma}
 
 \subsection{Plasma Core Power Balance}\label{plasma-core-power-balance}
 
@@ -1387,10 +1387,10 @@ The region directly outside the last closed flux surface of the core
 plasma is known as the scrape-off layer, and contains no structural
 material. Plasma entering this region is not confined and is removed by
 the divertor. PROCESS treats the scrape-off layer merely as a gap.
-Switch \texttt{iscrp} determines whether the inboard and outboard gaps
+Switch \texttt{i_plasma_wall_gap} determines whether the inboard and outboard gaps
 should be calculated as 10\% of the plasma minor radius
-(\texttt{iscrp\ =\ 0}), or set equal to the input values
-\texttt{scrapli} and \texttt{scraplo} (\texttt{iscrp\ =\ 1}).
+(\texttt{i_plasma_wall_gap\ =\ 0}), or set equal to the input values
+\texttt{scrapli} and \texttt{scraplo} (\texttt{i_plasma_wall_gap\ =\ 1}).
 
 \end{document}
 

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

@@ -78,17 +78,17 @@ derived directly from the energy confinement scaling law.
 
 `iradloss = 0` -- Total power lost is scaling power plus radiation:
 
-`pscaling + pradpv = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/plasma_volume`
+`pscaling + pradpv = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/vol_plasma`
 
 
 `iradloss = 1` -- Total power lost is scaling power plus radiation from a region defined as the "core":
 
-`pscaling + pcoreradpv = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/plasma_volume`
+`pscaling + pcoreradpv = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/vol_plasma`
 
 `iradloss = 2` -- Total power lost is scaling power only, with no additional 
 allowance for radiation. This is not recommended for power plant models.
 
-`pscaling = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/plasma_volume`
+`pscaling = f_alpha_plasma*alpha_power_density_total + charged_power_density + pden_plasma_ohmic_mw + pinjmw/vol_plasma`
 
 ## L-H Power Threshold Scalings
 

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

@@ -116,9 +116,9 @@ Switch value: `i_plasma_current = 1`
 
 The formula for calculating `fq` is:
 
-$$f_q = \left(\frac{{1.22 - 0.68  \epsilon}}{{(1.0 - \epsilon^2)^2}}\right)  \mathtt{{sf}}^2$$
+$$f_q = \left(\frac{{1.22 - 0.68  \epsilon}}{{(1.0 - \epsilon^2)^2}}\right)  \left(\frac{L_{\text{poloidal}}}{2\pi a}\right)^2$$
 
-Where $\epsilon$ is the inverse [aspect ratio](../plasma_geometry.md) ($\mathtt{eps}$) and $\mathtt{sf}$ is the shaping factor calculated in the [poloidal perimeter](../plasma_geometry.md#poloidal-perimeter) function in `plasma_geometry.py`
+Where $\epsilon$ is the inverse [aspect ratio](../plasma_geometry.md) ($\mathtt{eps}$) and $L_{\text{poloidal}}$ is the poloidal perimeter of the plasma.
 
 -----------
 
@@ -481,7 +481,7 @@ higher elongations causes an underestimate by up to 20%.
 
 Given the parameter dependencies a new plasma current relation based on fits to FIESTA
 equilibria was created. It showed no bias with any parameter fitted and that the fit is accurate to 10%.
-The linear relation between these and the 95% values expressed in does not hold at high values of elongation and triangularity as per the [`ishape = 0`](../plasma_geometry.md#plasma-geometry-parameters-geomty) relation.
+The linear relation between these and the 95% values expressed in does not hold at high values of elongation and triangularity as per the [`i_plasma_geometry = 0`](../plasma_geometry.md#plasma-geometry-parameters--plasma_geometry) relation.
 
 
 $$
@@ -539,10 +539,10 @@ $$
 For calculating the poloidal magnetic field created due to the presence of the plasma current, [Ampere's law](https://en.wikipedia.org/wiki/Amp%C3%A8re%27s_circuital_law) can be used. In this case the poloidal field is simply returned as:
 
 $$
-B_{\text{p}} = \frac{\mu_0 I_{\text{p}}}{\mathtt{pperim}}
+B_{\text{p}} = \frac{\mu_0 I_{\text{p}}}{\mathtt{len_plasma_poloidal}}
 $$
 
-Where `pperim` is the plasma poloidal perimeter calculated [here](../plasma_geometry.md#poloidal-perimeter).
+Where `len_plasma_poloidal` is the plasma poloidal perimeter calculated [here](../plasma_geometry.md#poloidal-perimeter).
 
 In the case of using the Peng double null scaling ([`i_plasma_current = 2`](plasma_current.md#star-peng-double-null-divertor-scaling-st)), the values $F_1$ and $F_2$ are calculated from [_plasc_bpol](plasma_current.md#_plasc_bpol) and used to calculated the poloidal field from the toroidal field as per:
 

documentation/proc-pages/scripts/plotting_scripts/plasma_cross_section_plot.py

@@ -0,0 +1,170 @@
+import math
+
+import numpy as np
+from bokeh.io import output_file, save
+from bokeh.layouts import column, row
+from bokeh.models import ColumnDataSource, CustomJS, Slider
+from bokeh.plotting import figure
+
+square = Slider(start=-1.0, end=1.0, value=0.0, step=0.01, title="Squareness")
+r0 = Slider(start=0.1, end=10, value=5.0, step=0.1, title="Major radius")
+a = Slider(start=0.01, end=10, value=3.0, step=0.01, title="Minor radius")
+delta = Slider(start=0.0, end=1, value=0.5, step=0.01, title="Triangularity | delta")
+kappa = Slider(start=0.0, end=10, value=2.0, step=0.01, title="Elongation | kappa")
+
+output_file("plass_shape_plot.html")
+x = np.linspace(-np.pi, np.pi, 256)
+
+# Sauter
+R = r0.value + a.value * np.cos(
+    x + delta.value * np.sin(x) - square.value * np.sin(2 * x)
+)
+Z = kappa.value * a.value * np.sin(x + square.value * np.sin(2 * x))
+
+# Mirror the Sauter plot
+R_mirror = -R
+
+# PROCESS
+x1 = (
+    2.0 * r0.value * (1.0 + delta.value)
+    - a.value * (delta.value**2 + kappa.value**2 - 1.0)
+) / (2.0 * (1.0 + delta.value))
+x2 = (
+    2.0 * r0.value * (delta.value - 1.0)
+    - a.value * (delta.value**2 + kappa.value**2 - 1.0)
+) / (2.0 * (delta.value - 1.0))
+r1 = 0.5 * math.sqrt(
+    (a.value**2 * ((delta.value + 1.0) ** 2 + kappa.value**2) ** 2)
+    / ((delta.value + 1.0) ** 2)
+)
+r2 = 0.5 * math.sqrt(
+    (a.value**2 * ((delta.value - 1.0) ** 2 + kappa.value**2) ** 2)
+    / ((delta.value - 1.0) ** 2)
+)
+theta1 = np.arcsin((kappa.value * a.value) / r1)
+theta2 = np.arcsin((kappa.value * a.value) / r2)
+
+inang = 1.0 / r1
+outang = 1.5 / r2
+
+angs1 = np.linspace(-theta1 + np.pi, (inang + theta1) + np.pi, 256, endpoint=True)
+angs2 = np.linspace(-(outang + theta2), theta2, 256, endpoint=True)
+
+xs1 = -(r1 * np.cos(angs1) - x1)
+ys1 = r1 * np.sin(angs1)
+xs2 = -(r2 * np.cos(angs2) - x2)
+ys2 = r2 * np.sin(angs2)
+
+# Mirror the PROCESS plot
+xs1_mirror = -xs1
+xs2_mirror = -xs2
+
+source = ColumnDataSource(
+    data={
+        "x": R,
+        "y": Z,
+        "x_mirror": R_mirror,
+        "y_mirror": Z,
+        "x1": xs1,
+        "y1": ys1,
+        "x2": xs2,
+        "y2": ys2,
+        "x1_mirror": xs1_mirror,
+        "y1_mirror": ys1,
+        "x2_mirror": xs2_mirror,
+        "y2_mirror": ys2,
+        "linspace": x,
+    }
+)
+
+plot = figure(
+    x_range=(-10, 10), y_range=(-10, 10), width=800, height=800, title="Plasma Shape"
+)
+plot.xaxis.axis_label = r"Radius"
+plot.yaxis.axis_label = r"Height"
+
+plot.line("x", "y", source=source, line_width=3, line_alpha=0.6, legend_label="Sauter")
+plot.line("x_mirror", "y_mirror", source=source, line_width=3, line_alpha=0.6)
+plot.line(
+    "x1",
+    "y1",
+    source=source,
+    line_width=3,
+    line_alpha=0.6,
+    color="red",
+    legend_label="PROCESS",
+)
+plot.line(
+    "x1_mirror", "y1_mirror", source=source, line_width=3, line_alpha=0.6, color="red"
+)
+plot.line("x2", "y2", source=source, line_width=3, line_alpha=0.6, color="red")
+plot.line(
+    "x2_mirror", "y2_mirror", source=source, line_width=3, line_alpha=0.6, color="red"
+)
+
+plot.legend.location = "top_left"
+
+callback = CustomJS(
+    args={
+        "source": source,
+        "r0": r0,
+        "a": a,
+        "delta": delta,
+        "kappa": kappa,
+        "square": square,
+    },
+    code="""
+    const A = r0.value;
+    const B = a.value;
+    const D = delta.value;
+    const K = kappa.value;
+    const S = square.value;
+    const x = source.data['linspace'];
+    const R = source.data['x'];
+    const Z = source.data['y'];
+    const R_mirror = source.data['x_mirror'];
+    const Z_mirror = source.data['y_mirror'];
+    const x1 = (2.0 * A * (1.0 + D) - B * (D**2 + K**2 - 1.0)) / (2.0 * (1.0 + D));
+    const x2 = (2.0 * A * (D - 1.0) - B * (D**2 + K**2 - 1.0)) / (2.0 * (D - 1.0));
+    const r1 = 0.5 * Math.sqrt((B**2 * ((D + 1.0) ** 2 + K**2) ** 2) / ((D + 1.0) ** 2));
+    const r2 = 0.5 * Math.sqrt((B**2 * ((D - 1.0) ** 2 + K**2) ** 2) / ((D - 1.0) ** 2));
+    const theta1 = Math.asin((K * B) / r1);
+    const theta2 = Math.asin((K * B) / r2);
+    const inang = 1.0 / r1;
+    const outang = 1.5 / r2;
+    const angs1 = Array.from({length: 256}, (_, i) => -theta1 + Math.PI + i * ((inang + 2*theta1) / 255));
+    const angs2 = Array.from({length: 256}, (_, i) => -(outang + theta2) + i * ((theta2 + (outang + theta2)) / 255));
+    const xs1 = source.data['x1'];
+    const ys1 = source.data['y1'];
+    const xs2 = source.data['x2'];
+    const ys2 = source.data['y2'];
+    const xs1_mirror = source.data['x1_mirror'];
+    const ys1_mirror = source.data['y1_mirror'];
+    const xs2_mirror = source.data['x2_mirror'];
+    const ys2_mirror = source.data['y2_mirror'];
+    for (let i = 0; i < x.length; i++) {
+        R[i] = A + B * Math.cos(x[i] + D * Math.sin(x[i]) - S * Math.sin(2 * x[i]));
+        Z[i] = K * B * Math.sin(x[i] + S * Math.sin(2 * x[i]));
+        R_mirror[i] = -R[i];
+        Z_mirror[i] = Z[i];
+        xs1[i] = -(r1 * Math.cos(angs1[i]) - x1);
+        ys1[i] = r1 * Math.sin(angs1[i]);
+        xs2[i] = -(r2 * Math.cos(angs2[i]) - x2);
+        ys2[i] = r2 * Math.sin(angs2[i]);
+        xs1_mirror[i] = -xs1[i];
+        ys1_mirror[i] = ys1[i];
+        xs2_mirror[i] = -xs2[i];
+        ys2_mirror[i] = ys2[i];
+    }
+    source.change.emit();
+""",
+)
+
+r0.js_on_change("value", callback)
+a.js_on_change("value", callback)
+delta.js_on_change("value", callback)
+kappa.js_on_change("value", callback)
+square.js_on_change("value", callback)
+
+# Save the plot as HTML
+save(row(plot, column(r0, a, delta, kappa, square)), filename="./plasma_plot.html")