Remove gtscale and enhance iprofile

documentation/proc-pages/index.md

@@ -7,7 +7,7 @@
     class="logo"
     >
 
-PROCESS is a systems code at [CCFE](https://ccfe.ukaea.uk/) that calculates in a
+PROCESS is a systems code at [UKAEA](https://ccfe.ukaea.uk/) that calculates in a
 self-consistent manner the parameters of a fusion power plant with a specified 
 performance, ensuring that its operating limits are not violated, and with the option 
 to optimise to a given function of these parameters.
@@ -72,8 +72,6 @@ to anyone using PROCESS outputs or models based on them.
 - [Stuart Muldrew](mailto:stuart.muldrew@ukaea.uk)
 - [Alex Pearce](mailto:alexander.pearce@ukaea.uk)
 - [Jonathan Maddock](mailto:jonathan.maddock@ukaea.uk)
-- [Charles Griesel](mailto:charles.griesel@ukaea.uk)
-- [Rhian Chapman](mailto:rhian.chapman@ukaea.uk)
 - [Timothy Nunn](mailto:timothy.nunn@ukaea.uk)
 - [Christopher Ashe](mailto:christopher.ashe@ukaea.uk)
 - [Georgina Graham](mailto:georgina.graham@ukaea.uk)

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

@@ -3,7 +3,7 @@
 The plasma beta limit[^1] is given by 
 
 $$\begin{aligned}
-\beta < g \, \frac{I(\mbox{MA})}{a(\mbox{m}) \, B_0(\mbox{T})}
+\beta < 0.01\, g \, \frac{I(\mbox{MA})}{a(\mbox{m}) \, B_0(\mbox{T})}
 \end{aligned}$$
 
 where $B_0$ is the axial vacuum toroidal field. The beta
@@ -20,19 +20,21 @@ By default, $\beta$ is defined with respect to the total equilibrium B-field [^2
 | 2 | Apply the $\beta$ limit to only the thermal plus neutral beam contributions to beta |
 | 3 | Apply the $\beta$ limit to the total beta (including the contribution from fast ions), calculated using only the toroidal field |
 
-### Scaling of beta $g$ coefficient
+### Setting the Beta $g$ Coefficient
 
-Switch `gtscale` determines how the beta $g$ coefficient `dnbeta` should 
-be calculated, using the inverse aspect ratio $\epsilon = a/R$.
+Switch `iprofile` determines how the beta $g$ coefficient `dnbeta` should 
+be calculated.
 
-| `gtscale` | Description |
+| `iprofile` | Description |
 | :-: | - |
-| 0 | `dnbeta` is an input. |
-| 1 | $g=2.7(1+5\epsilon^{3.5})$ (which gives g = 3.0 for aspect ratio = 3) |
-| 2 | $g=3.12+3.5\epsilon^{1.7}$ (based on Menard et al. "Fusion Nuclear Science Facilities and Pilot Plants Based on the Spherical Tokamak", Nucl. Fusion, 2016, 44)  |
+| 0 | `alphaj`, `rli` and `dnbeta` are inputs. |
+| 1 | `alphaj`, `rli` and `dnbeta` are calulcated consistently. `dnbeta` calculated using $g=4l_i$ [^3].  This is only recommended for high aspect ratio tokamaks.|
+| 2 | `alphaj` and `rli` are inputs. `dnbeta` calculated using $g=2.7(1+5\epsilon^{3.5})$ (which gives g = 3.0 for aspect ratio = 3) |
+| 3 | `alphaj` and `rli` are inputs. `dnbeta` calculated using $g=3.12+3.5\epsilon^{1.7}$ [^4]|
+| 4 | `alphaj` and `dnbeta` are inputs. `rli` calculated from elongation [^4]. This is only recommended for spherical tokamaks.|
+| 5 | `alphaj` is an input.  `rli` calculated from elongation and `dnbeta` calculated using $g=3.12+3.5\epsilon^{1.7}$ [^4]. This is only recommended for spherical tokamaks.|
 
-!!! Note 
-    `gtscale` is over-ridden if `iprofile` = 1.
+Further details on the calculation of `alphaj` and `rli` is given in [Plasma Current](./plasma_current.md).
 
 ### Limiting $\epsilon\beta_p$
 
@@ -43,4 +45,8 @@ is be set using input parameter `epbetmax`.
 
 [^1]: N.A. Uckan and ITER Physics Group, 'ITER Physics Design Guidelines: 1989',
 
-[^2]: D.J. Ward, 'PROCESS Fast Alpha Pressure', Work File Note F/PL/PJK/PROCESS/CODE/050
+[^2]: D.J. Ward, 'PROCESS Fast Alpha Pressure', Work File Note F/PL/PJK/PROCESS/CODE/050
+
+[^3]: Tokamaks 4th Edition, Wesson, page 116
+
+[^4]: Menard et al. (2016), Nuclear Fusion, 56, 106023

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

@@ -1,6 +1,6 @@
 # Plasma Current Scaling Laws
 
-A number of plasma current scaling laws are available in PROCESS $[^1]. These are calculated in 
+A number of plasma current scaling laws are available in PROCESS [^1]. These are calculated in 
 routine `culcur`, which is called by `physics`. The safety factor $q_{95}$ required to prevent 
 disruptive MHD instabilities dictates the plasma current Ip:
 
@@ -28,21 +28,31 @@ the switch `icurr`, as follows:
 
 A limited degree of self-consistency between the plasma current profile and other parameters [^6] can be 
 enforced by setting switch `iprofile = 1`. This sets the current 
-profile peaking factor $\alpha_J$ (`alphaj`) and the normalised internal inductance $l_i$ (`rli`) using the 
+profile peaking factor $\alpha_J$ (`alphaj`),  the normalised internal inductance $l_i$ (`rli`) and beta limit $g$-factor (`dnbeta`) using the 
 safety factor on axis `q0` and the cylindrical safety factor $q*$ (`qstar`):   
 
 $$\begin{aligned}
 \alpha_J = \frac{q*}{q_0} - 1
 \end{aligned}$$
 
 $$\begin{aligned}
-l_i = ln(1.65+0.89\alpha_J)
+l_i = \rm{ln}(1.65+0.89\alpha_J)
 \end{aligned}$$
 
-The beta $g$ coefficient `dnbeta` also scales with $l_i$, as described above.
+$$\begin{aligned}
+g = 4 l_i
+\end{aligned}$$
 
 It is recommended that current scaling law `icurr = 4` is used if `iprofile = 1`. 
-Switch `gtscale` is over-ridden if `iprofile = 1`.
+This relation is only applicable to large aspect ratio tokamaks.
+
+For spherical tokamaks, the normalised internal inductance can be set from the elongation using `iprofile = 4` or `iprofile = 5`:
+
+$$\begin{aligned}
+l_i = 3.4 - \kappa_x
+\end{aligned}$$
+
+Further desciption of `iprofile` is given in [Beta Limit](./plasma_beta.md).
 
 ## Bootstrap, Diamagnetic and Pfirsch-Schlüter Current Scalings
 
@@ -56,8 +66,8 @@ existence of pedestals, whereas the Sauter et al. scaling
 | :-: | - |
 | 1 | ITER scaling -- To use the ITER scaling method for the bootstrap current fraction.  Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$). This method is valid at high aspect ratio only.
 | 2 | General scaling -- To use a more general scaling method, set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
-| 3 | Numerically fitted scaling [^7] -- To use a numerically fitted scaling method, valid for all aspect ratios, set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
-| 4 | Sauter, Angioni and Lin-Liu scaling [^8] [^9] -- Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
+| 3 | Numerically fitted scaling [^8] -- To use a numerically fitted scaling method, valid for all aspect ratios, set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
+| 4 | Sauter, Angioni and Lin-Liu scaling [^9] [^10] -- Set `bscfmax` to the maximum required bootstrap current fraction ($\leq 1$).
 
 !!! Note "Fixed Bootstrap Current"
     Direct input -- To input the bootstrap current fraction directly, set `bscfmax` 
@@ -101,7 +111,7 @@ Unpublished internal Oak Ridge document.
 Current Drive', ITER-TN-PH-8-4, 13--17 June 1988, Garching, FRG
 [^6]: Y. Sakamoto, 'Recent progress in vertical stability analysis in JA',
 Task meeting EU-JA #16, Fusion for Energy, Garching, 24--25 June 2014
-
-[^7]: H.R. Wilson, Nuclear Fusion **32** (1992) 257
-[^8]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **6** (1999) 2834 
-[^9]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **9** (2002) 5140
+[^7]: Menard et al. (2016), Nuclear Fusion, 56, 106023
+[^8]: H.R. Wilson, Nuclear Fusion **32** (1992) 257
+[^9]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **6** (1999) 2834 
+[^10]: O. Sauter, C. Angioni and Y.R. Lin-Liu, Physics of Plasmas **9** (2002) 5140

process/physics.py

@@ -1381,28 +1381,26 @@ def physics(self):
         )
 
         # Calculate physics_variables.beta limit
-        if physics_variables.iprofile == 0:
-            if physics_variables.gtscale == 1:
-                # Original scaling law
-                physics_variables.dnbeta = 2.7e0 * (
-                    1.0e0 + 5.0e0 * physics_variables.eps**3.5e0
-                )
-
-            if physics_variables.gtscale == 2:
-                # See Issue #1439
-                # physics_variables.dnbeta found from physics_variables.aspect ratio scaling on p32 of Menard:
-                # Menard, et al. "Fusion Nuclear Science Facilities
-                # and Pilot Plants Based on the Spherical Tokamak."
-                # Nucl. Fusion, 2016, 44.
-                physics_variables.dnbeta = (
-                    3.12e0 + 3.5e0 * physics_variables.eps**1.7e0
-                )
 
-        else:
+        if physics_variables.iprofile == 1:
             # Relation between physics_variables.beta limit and plasma internal inductance
             # Hartmann and Zohm
             physics_variables.dnbeta = 4.0e0 * physics_variables.rli
 
+        if physics_variables.iprofile == 2:
+            # Original scaling law
+            physics_variables.dnbeta = 2.7e0 * (
+                1.0e0 + 5.0e0 * physics_variables.eps**3.5e0
+            )
+
+        if physics_variables.iprofile == 3 or physics_variables.iprofile == 5:
+            # physics_variables.dnbeta found from physics_variables.aspect ratio scaling on p32 of Menard:
+            # Menard, et al. "Fusion Nuclear Science Facilities
+            # and Pilot Plants Based on the Spherical Tokamak."
+            # Nucl. Fusion, 2016, 44.
+            physics_variables.dnbeta = 3.12e0 + 3.5e0 * physics_variables.eps**1.7e0
+
+        # culblm returns the betalim for beta
         physics_variables.betalim = culblm(
             physics_variables.bt,
             physics_variables.dnbeta,
@@ -2030,8 +2028,12 @@ def culcur(
         8 = Sauter scaling (allowing negative triangularity) Issue #392
             'Geometric formulas for system codes including the effect of negative triangularity'
         iprofile : input integer : switch for current profile consistency
-        0 = use input alphaj, rli
-        1 = make these consistent with q, q0
+        0 use input values for alphaj, rli, dnbeta
+        1 make these consistent with input q, q_0 values (recommend `icurr=4` with this option)
+        2 use input values for alphaj, rli. Scale dnbeta with aspect ratio (original scaling)
+        3 use input values for alphaj, rli. Scale dnbeta with aspect ratio (Menard scaling)
+        4 use input values for alphaj, dnbeta. Set rli from elongation (Menard scaling)
+        5 use input value for alphaj.  Set rli and dnbeta from Menard scaling
         kappa    : input real :  plasma elongation
         kappa95  : input real :  plasma elongation at 95% surface
         p0       : input real :  central plasma pressure (Pa)
@@ -2151,12 +2153,18 @@ def culcur(
             icurr, plascur, q95, asp, eps, bt, kappa, triang, pperim, constants.rmu0
         )
 
-        # Ensure current profile consistency, if required
-        # This is as described in Hartmann and Zohm only if icurr = 4 as well...
-
         if iprofile == 1:
+            # Ensure current profile consistency, if required
+            # This is as described in Hartmann and Zohm only if icurr = 4 as well...
+
+            # Tokamaks 4th Edition, Wesson, page 116
             alphaj = qstar / q0 - 1.0
-            rli = np.log(1.65 + 0.89 * alphaj)  # Tokamaks 4th Edition, Wesson, page 116
+            rli = np.log(1.65 + 0.89 * alphaj)
+
+        if iprofile == 4 or iprofile == 5:
+            # Spherical Tokamak relation for internal inductance
+            # Menard et al. (2016), Nuclear Fusion, 56, 106023
+            rli = 3.4 - kappa
 
         return alphaj, rli, bp, qstar, plascur
 
@@ -2494,15 +2502,15 @@ def outplas(self):
             po.osubhd(self.outfile, "Current and Field :")
 
             if stellarator_variables.istell == 0:
-                if physics_variables.iprofile == 0:
+                if physics_variables.iprofile == 1:
                     po.ocmmnt(
                         self.outfile,
-                        "Consistency between q0,q,alphaj,rli,dnbeta is not enforced",
+                        "Consistency between q0,q,alphaj,rli,dnbeta is enforced",
                     )
                 else:
                     po.ocmmnt(
                         self.outfile,
-                        "Consistency between q0,q,alphaj,rli,dnbeta is enforced",
+                        "Consistency between q0,q,alphaj,rli,dnbeta is not enforced",
                     )
 
                 po.oblnkl(self.outfile)

source/fortran/input.f90

@@ -309,7 +309,7 @@ subroutine parse_input_file(in_file,out_file,show_changes)
       fgwsep, rhopedn, tratio, q0, ishape, fne0, ignite, ftrit, &
       ifalphap, tauee_in, alphaj, alphat, icurr, q, ti, tesep, rli, triang, &
       itart, ralpne, iprofile, triang95, rad_fraction_sol, betbm0, protium, &
-      teped, fhe3, iwalld, gamma, falpha, fgwped, gtscale, tbeta, ibss, &
+      teped, fhe3, iwalld, gamma, falpha, fgwped, tbeta, ibss, &
       iradloss, te, alphan, rmajor, kappa, iinvqd, fkzohm, beamfus0, &
       tauratio, idensl, ieped, bt, iscrp, ipnlaws, betalim, betalim_lower, &
       idia, ips, m_s_limit, burnup_in
@@ -612,9 +612,6 @@ subroutine parse_input_file(in_file,out_file,show_changes)
        case ('gamma')
           call parse_real_variable('gamma', gamma, 0.1D0, 1.0D0, &
                'Ejima coefficient for resistive V-s formula')
-       case ('gtscale')
-          call parse_int_variable('gtscale', gtscale, 0, 2, &
-               'Flag to scale beta coefficient with R/a')
        case ('hfact')
           call parse_real_variable('hfact', hfact, 0.01D0, 10.0D0, &
                'Energy confinement time H factor')
@@ -653,7 +650,7 @@ subroutine parse_input_file(in_file,out_file,show_changes)
           call parse_int_variable('ipedestal', ipedestal, 0, 1, &
                'Switch for plasma profile type')
        case ('iprofile')
-          call parse_int_variable('iprofile', iprofile, 0, 1, &
+          call parse_int_variable('iprofile', iprofile, 0, 5, &
                'Switch for current profile consistency')
        case ('ips')
           call parse_int_variable('ips', ips, 0, 1, &

source/fortran/physics_variables.f90

@@ -27,7 +27,7 @@ module physics_variables
   !! average mass of all ions (amu)
 
   real(dp) :: alphaj
-  !! current profile index (calculated from q_0, q if `iprofile=1`)
+  !! current profile index (calculated from q_0 and q if `iprofile=1`)
 
   real(dp) :: alphan
   !! density profile index
@@ -124,8 +124,7 @@ module physics_variables
   !! hot beam ion density from calculation (/m3)
 
   real(dp) :: dnbeta
-  !! Troyon-like coefficient for beta scaling calculated
-  !! as 4*rli if `iprofile=1` (see also gtscale option)
+  !! Troyon-like coefficient for beta scaling
 
   real(dp) :: dnelimt
   !! density limit (/m3)
@@ -231,13 +230,6 @@ module physics_variables
   real(dp) :: gammaft
   !! ratio of (fast alpha + neutral beam beta) to thermal beta
 
-  integer :: gtscale
-  !! switch for a/R scaling of dnbeta (`iprofile=0` only):
-  !!
-  !! - =0 do not scale dnbeta with eps
-  !! - =1 scale dnbeta with eps, original scaling
-  !! - =2 scale dnbeta with eps, Menard scaling
-
   real(dp), dimension(ipnlaws) :: hfac
   !! H factors for an ignited plasma for each energy confinement time scaling law
 
@@ -379,8 +371,12 @@ module physics_variables
   integer :: iprofile
   !! switch for current profile consistency:
   !!
-  !! - =0 use input values for alphaj, rli, dnbeta (but see gtscale option)
+  !! - =0 use input values for alphaj, rli, dnbeta
   !! - =1 make these consistent with input q, q_0 values (recommend `icurr=4` with this option)
+  !! - =2 use input values for alphaj, rli. Scale dnbeta with aspect ratio (original scaling)
+  !! - =3 use input values for alphaj, rli. Scale dnbeta with aspect ratio (Menard scaling)
+  !! - =4 use input values for alphaj, dnbeta. Set rli from elongation (Menard scaling)
+  !! - =5 use input value for alphaj.  Set rli and dnbeta from Menard scaling
 
   integer :: iradloss
   !! switch for radiation loss term usage in power balance (see User Guide):
@@ -954,7 +950,6 @@ subroutine init_physics_variables
     fvsbrnni = 1.0D0
     gamma = 0.4D0
     gammaft = 0.0D0
-    gtscale = 0
     hfac = 0.0D0
     hfact = 1.0D0
     taumax = 10.0D0