From a6c665223469b4e83fc8cac9535701488082a1f6 Mon Sep 17 00:00:00 2001 From: james <81617086+je-cook@users.noreply.github.com> Date: Mon, 4 Mar 2024 07:32:30 +0000 Subject: [PATCH] Numba-ify physics functions --- process/physics.py | 5066 ++++++++++++++++++------------------ tests/unit/test_physics.py | 54 +- 2 files changed, 2552 insertions(+), 2568 deletions(-) diff --git a/process/physics.py b/process/physics.py index 0bfc45a4f6..b09fff0af0 100644 --- a/process/physics.py +++ b/process/physics.py @@ -1,4 +1,5 @@ import numpy as np +import numba as nb import scipy.integrate as integrate from scipy.optimize import root_scalar import math @@ -28,6 +29,753 @@ ) +@nb.jit(nopython=True, cache=True) +def coulg(j, tempe, ne): + """Coulomb logarithm + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + This function calculates the Coulomb logarithm, valid + for e-e collisions (T_e > 0.01 keV), and for + e-i collisions (T_e > 0.01*Zeff^2) (Alexander, 9/5/1994). +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + C. A. Ordonez and M. I. Molina, Phys. Plasmas 1 (1994) 2515 + Rev. Mod. Phys., V.48, Part 1 (1976) 275 + """ + return 15.9 - 0.5 * np.log(ne[j - 1]) + np.log(tempe[j - 1]) + + +@nb.jit(nopython=True, cache=True) +def nuee(j, tempe, ne): + """Frequency of electron-electron collisions + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + This function calculates the frequency of electron-electron + collisions (Hz): NUEE = 4*SQRT(pi)/3*Ne*e**4*lambd/ + SQRT(Me)/Te**1.5 +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + Yushmanov, 25th April 1987 (?), + updated by Pereverzev, 9th November 1994 (?) + """ + return ( + 670.0 * coulg(j, tempe, ne) * ne[j - 1] / (tempe[j - 1] * np.sqrt(tempe[j - 1])) + ) + + +@nb.jit(nopython=True, cache=True) +def nui(j, zmain, ni, tempi, amain): + """Full frequency of ion collisions + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + This function calculates the full frequency of ion + collisions (Hz). +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + None + """ + # Coulomb logarithm = 15 is used + return ( + zmain[j - 1] ** 4 + * ni[j - 1] + * 322.0 + / (tempi[j - 1] * np.sqrt(tempi[j - 1] * amain[j - 1])) + ) + + +@nb.jit(nopython=True, cache=True) +def _plascar_bpol(aspect, kappa, delta): + # Original coding, only suitable for TARTs [STAR Code] + + eps = 1.0 / aspect + + c1 = kappa**2 / (1.0 + delta) + delta + c2 = kappa**2 / (1.0 - delta) - delta + + d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 + d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 + + c1_aspect = (c1 * eps - 1.0) if aspect < c1 else (1.0 - c1 * eps) + + y1 = np.sqrt(c1_aspect / (1.0 + eps)) * (1.0 + delta) / kappa + y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * (1.0 - delta) / kappa + + h2 = (1.0 + (c2 - 1.0) * eps / 2.0) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) + f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * (c2 * eps + 1.0)) + g = eps * kappa / (1.0 - eps * delta) + ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) + + h1 = (1.0 + (1.0 - c1) * eps / 2.0) / np.sqrt((1.0 + eps) * c1_aspect) + f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) + + if aspect < c1: + ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) + else: + ff1 = -f1 * (-g + 2.0 * h1 * np.arctan(y1)) + + return ff1, ff2, d1, d2 + + +@nb.jit(nopython=True, cache=True) +def bpol(icurr, ip, qbar, aspect, bt, kappa, delta, perim, rmu0): + """Function to calculate poloidal field + author: J Galambos, FEDC/ORNL + author: P J Knight, CCFE, Culham Science Centre + icurr : input integer : current scaling model to use + ip : input real : plasma current (A) + qbar : input real : edge q-bar + aspect : input real : plasma aspect ratio + bt : input real : toroidal field on axis (T) + kappa : input real : plasma elongation + delta : input real : plasma triangularity + perim : input real : plasma perimeter (m) + This function calculates the poloidal field in Tesla, + using a simple calculation using Stoke's Law for conventional + tokamaks, or for TARTs, a scaling from Peng, Galambos and + Shipe (1992). + J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, + Fusion Technology, 21, 1729 + """ + if icurr != 2: + return rmu0 * ip / perim + + ff1, ff2, _, _ = _plascar_bpol(aspect, kappa, delta) + + return bt * (ff1 + ff2) / (2.0 * np.pi * qbar) + + +@nb.jit(nopython=True, cache=True) +def plasc(qbar, aspect, rminor, bt, kappa, delta): + """Function to calculate plasma current (Peng scaling) + author: J Galambos, FEDC/ORNL + author: P J Knight, CCFE, Culham Science Centre + aspect : input real : plasma aspect ratio + bt : input real : toroidal field on axis (T) + delta : input real : plasma triangularity + kappa : input real : plasma elongation + qbar : input real : edge q-bar + rminor : input real : plasma minor radius (m) + This function calculates the plasma current in MA, + using a scaling from Peng, Galambos and Shipe (1992). + It is primarily used for Tight Aspect Ratio Tokamaks and is + selected via icurr=2. + J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, + unpublished internal Oak Ridge document + Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, + Fusion Technology, 21, 1729 + """ + + ff1, ff2, d1, d2 = _plascar_bpol(aspect, kappa, delta) + + e1 = 2.0 * kappa / (d1 * (1.0 + delta)) + e2 = 2.0 * kappa / (d2 * (1.0 - delta)) + + return ( + rminor + * bt + / qbar + * 5.0 + * kappa + / (2.0 * np.pi**2) + * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) + * (ff1 + ff2) + ) + + +@nb.jit(nopython=True, cache=True) +def rether(alphan, alphat, dene, dlamie, te, ti, zeffai): + """Routine to find the equilibration power between the + ions and electrons + author: P J Knight, CCFE, Culham Science Centre + alphan : input real : density profile index + alphat : input real : temperature profile index + dene : input real : electron density (/m3) + dlamie : input real : ion-electron coulomb logarithm + te : input real : electron temperature (keV) + ti : input real : ion temperature (keV) + zeffai : input real : mass weighted plasma effective charge + piepv : output real : ion/electron equilibration power (MW/m3) + This routine calculates the equilibration power between the + ions and electrons. + Unknown origin + """ + profie = (1.0 + alphan) ** 2 / ( + (2.0 * alphan - 0.5 * alphat + 1.0) * np.sqrt(1.0 + alphat) + ) + conie = 2.42165e-41 * dlamie * dene**2 * zeffai * profie + + return conie * (ti - te) / (te**1.5) + + +@nb.jit(nopython=True, cache=True) +def vscalc( + csawth, + eps, + facoh, + gamma, + kappa, + rmajor, + rplas, + plascur, + t_fusion_ramp, + tburn, + rli, + rmu0, +): + """Volt-second requirements + author: P J Knight, CCFE, Culham Science Centre + csawth : input real : coefficient for sawteeth effects + eps : input real : inverse aspect ratio + facoh : input real : fraction of plasma current produced inductively + gamma : input real : Ejima coeff for resistive start-up V-s component + kappa : input real : plasma elongation + plascur: input real : plasma current (A) + rli : input real : plasma normalised inductivity + rmajor : input real : plasma major radius (m) + rplas : input real : plasma resistance (ohm) + t_fusion_ramp : input real : heating time (s) + tburn : input real : burn time (s) + phiint : output real : internal plasma volt-seconds (Wb) + rlp : output real : plasma inductance (H) + vsbrn : output real : volt-seconds needed during flat-top (heat+burn) (Wb) + vsind : output real : internal and external plasma inductance V-s (Wb) + vsres : output real : resistive losses in start-up volt-seconds (Wb) + vsstt : output real : total volt-seconds needed (Wb) + This subroutine calculates the volt-second requirements and some + other related items. + AEA FUS 251: A User's Guide to the PROCESS Systems Code + """ + # Internal inductance + + rlpint = rmu0 * rmajor * rli / 2.0 + phiint = rlpint * plascur + + # Start-up resistive component + # Uses ITER formula without the 10 V-s add-on + + vsres = gamma * rmu0 * plascur * rmajor + + # Hirshman, Neilson: Physics of Fluids, 29 (1986) p790 + # fit for external inductance + + aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( + 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps + ) + beps = ( + 0.73 * np.sqrt(eps) * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) + ) + rlpext = rmajor * rmu0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) + + rlp = rlpext + rlpint + + # Inductive V-s component + + vsind = rlp * plascur + vsstt = vsres + vsind + + # Loop voltage during flat-top + # Include enhancement factor in flattop V-s requirement + # to account for MHD sawtooth effects. + + vburn = plascur * rplas * facoh * csawth + + # N.B. tburn on first iteration will not be correct + # if the pulsed reactor option is used, but the value + # will be correct on subsequent calls. + + vsbrn = vburn * (t_fusion_ramp + tburn) + vsstt = vsstt + vsbrn + + return phiint, rlp, vsbrn, vsind, vsres, vsstt + + +@nb.jit(nopython=True, cache=True) +def culblm(bt, dnbeta, plascur, rminor): + """Beta scaling limit + author: P J Knight, CCFE, Culham Science Centre + bt : input real : toroidal B-field on plasma axis (T) + dnbeta : input real : Troyon-like g coefficient + plascur : input real : plasma current (A) + rminor : input real : plasma minor axis (m) + betalim : output real : beta limit as defined below + This subroutine calculates the beta limit, using + the algorithm documented in AEA FUS 172. + The limit applies to beta defined with respect to the total B-field. + Switch iculbl determines which components of beta to include. + + If iculbl = 0, then the limit is applied to the total beta + If iculbl = 1, then the limit is applied to the thermal beta only + If iculbl = 2, then the limit is applied to the thermal + neutral beam beta components + If iculbl = 3, then the limit is applied to the toroidal beta + + The default value for the g coefficient is dnbeta = 3.5 + AEA FUS 172: Physics Assessment for the European Reactor Study + AEA FUS 251: A User's Guide to the PROCESS Systems Code + """ + + return 0.01 * dnbeta * (plascur / 1.0e6) / (rminor * bt) + + +@nb.jit(nopython=True, cache=True) +def conhas(alphaj, alphap, bt, delta95, eps, kappa95, p0, rmu0): + """Routine to calculate the F coefficient used for scaling the + plasma current + author: P J Knight, CCFE, Culham Science Centre + alphaj : input real : current profile index + alphap : input real : pressure profile index + bt : input real : toroidal field on axis (T) + delta95 : input real : plasma triangularity 95% + eps : input real : inverse aspect ratio + kappa95 : input real : plasma elongation 95% + p0 : input real : central plasma pressure (Pa) + fq : output real : scaling for edge q from circular + cross-section cylindrical case + This routine calculates the F coefficient used for scaling the + plasma current, using the Connor-Hastie scaling given in + AEA FUS 172. + AEA FUS 172: Physics Assessment for the European Reactor Study + AEA FUS 251: A User's Guide to the PROCESS Systems Code + """ + + # Exponent in Connor-Hastie current profile - matching total + # current gives the following trivial relation + lamda = alphaj + + # Exponent in Connor-Hastie pressure profile + nu = alphap + + # Central plasma beta + beta0 = 2.0 * rmu0 * p0 / (bt**2) + + # Plasma internal inductance + lamp1 = 1.0 + lamda + li = lamp1 / lamda * (lamp1 / lamda * np.log(lamp1) - 1.0) + + # T/r in AEA FUS 172 + kap1 = kappa95 + 1.0 + tr = kappa95 * delta95 / kap1**2 + + # E/r in AEA FUS 172 + er = (kappa95 - 1.0) / kap1 + + # T primed in AEA FUS 172 + tprime = 2.0 * tr * lamp1 / (1.0 + 0.5 * lamda) + + # E primed in AEA FUS 172 + eprime = er * lamp1 / (1.0 + lamda / 3.0) + + # Delta primed in AEA FUS 172 + deltap = (0.5 * kap1 * eps * 0.5 * li) + ( + beta0 / (0.5 * kap1 * eps) + ) * lamp1**2 / (1.0 + nu) + + # Delta/R0 in AEA FUS 172 + deltar = beta0 / 6.0 * (1.0 + 5.0 * lamda / 6.0 + 0.25 * lamda**2) + ( + 0.5 * kap1 * eps + ) ** 2 * 0.125 * (1.0 - (lamda**2) / 3.0) + + # F coefficient + return (0.5 * kap1) ** 2 * ( + 1.0 + + eps**2 * (0.5 * kap1) ** 2 + + 0.5 * deltap**2 + + 2.0 * deltar + + 0.5 * (eprime**2 + er**2) + + 0.5 * (tprime**2 + 4.0 * tr**2) + ) + + +def bsinteg( + y, + dene, + ten, + bt, + rminor, + rmajor, + zeff, + alphat, + alphan, + q0, + q95, + betat, + echarge, + rmu0, +): + """Integrand function for Nevins et al bootstrap current scaling + author: P J Knight, CCFE, Culham Science Centre + y : input real : abscissa of integration, = normalised minor radius + This function calculates the integrand function for the + Nevins et al bootstrap current scaling, 4/11/90. + AEA FUS 251: A User's Guide to the PROCESS Systems Code + """ + # Constants for fit to q-profile + c1 = 1.0 + c2 = 1.0 + c3 = 1.0 + + # Compute average electron beta + betae = dene * ten * 1.0e3 * echarge / (bt**2 / (2.0 * rmu0)) + + nabla = rminor * np.sqrt(y) / rmajor + x = (1.46 * np.sqrt(nabla) + 2.4 * nabla) / (1.0 - nabla) ** 1.5 + z = zeff + d = ( + 1.414 * z + + z * z + + x * (0.754 + 2.657 * z + 2.0 * z * z) + + x * x * (0.348 + 1.243 * z + z * z) + ) + al2 = -x * (0.884 + 2.074 * z) / d + a2 = alphat * (1.0 - y) ** (alphan + alphat - 1.0) + alphai = -1.172 / (1.0 + 0.462 * x) + a1 = (alphan + alphat) * (1.0 - y) ** (alphan + alphat - 1.0) + al1 = x * (0.754 + 2.21 * z + z * z + x * (0.348 + 1.243 * z + z * z)) / d + + # q-profile + q = q0 + (q95 - q0) * (c1 * y + c2 * y * y + c3 * y**3) / (c1 + c2 + c3) + + pratio = (betat - betae) / betae + + return (q / q95) * (al1 * (a1 + pratio * (a1 + alphai * a2)) + al2 * a2) + + +def diamagnetic_fraction_hender(beta): + """author: S.I. Muldrew, CCFE, Culham Science Centre + Diamagnetic contribution at tight aspect ratio. + Tim Hender fit + """ + return beta / 2.8 + + +def diamagnetic_fraction_scene(beta, q95, q0): + """author: S.I. Muldrew, CCFE, Culham Science Centre + Diamagnetic fraction based on SCENE fit by Tim Hender + See Issue #992 + """ + return beta * (0.1 * q95 / q0 + 0.44) * 4.14e-1 + + +def ps_fraction_scene(beta): + """author: S.I. Muldrew, CCFE, Culham Science Centre + Pfirsch-Schlüter fraction based on SCENE fit by Tim Hender + See Issue #992 + """ + return -9e-2 * beta + + +@nb.jit(nopython=True, cache=True) +def nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni): + """Relative frequency of ion collisions + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + This function calculates the relative frequency of ion + collisions: NU* = Nui*q*Rt/eps**1.5/Vti + The full ion collision frequency NUI is used. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + Yushmanov, 30th April 1987 (?) + """ + return ( + 3.2e-6 + * nui(j, zmain, ni, tempi, amain) + * rmajor + / ( + np.abs(mu[j - 1] + 1.0e-4) + * sqeps[j - 1] ** 3 + * np.sqrt(tempi[j - 1] / amain[j - 1]) + ) + ) + + +@nb.jit(nopython=True, cache=True) +def dcsa(j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps): + """Grad(ln(ne)) coefficient in the Sauter bootstrap scaling + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 2 to nr + nr : input integer : maximum value of j + This function calculates the coefficient scaling grad(ln(ne)) + in the Sauter bootstrap current scaling. + Code by Angioni, 29th May 2002. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + O. Sauter, C. Angioni and Y. R. Lin-Liu, + Physics of Plasmas 6 (1999) 2834 + O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) + Physics of Plasmas 9 (2002) 5140 + + DCSA $\\equiv \\mathcal{L}_{31}$, Eq.14a, Sauter et al, 1999 + """ + zz = zef[j - 1] + zft = tpf(j, triang, sqeps) + _nues = nues(j, rmajor, zef, mu, sqeps, tempe, ne) + zdf = (1.0 + (1.0 - 0.1 * zft) * np.sqrt(_nues)) + (0.5 * (1.0 - zft) * _nues) / zz + zft /= zdf # $f^{31}_{teff}(\nu_{e*})$, Eq.14b + zft2 = zft**2 + zft3 = zft2 * zft + dcsa = ( + ((1.0 + 1.4 / (zz + 1.0)) * zft) + - (1.9 / (zz + 1.0) * zft2) + + ((0.3 * zft3 + 0.2 * zft3 * zft) / (zz + 1.0)) + ) + + # Corrections suggested by Fable, 15/05/2015 + return dcsa * beta_poloidal_local_total( + j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho + ) + + +@nb.jit(nopython=True, cache=True) +def hcsa(j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps): + """Grad(ln(Te)) coefficient in the Sauter bootstrap scaling + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 2 to nr + nr : input integer : maximum value of j + This function calculates the coefficient scaling grad(ln(Te)) + in the Sauter bootstrap current scaling. + Code by Angioni, 29th May 2002. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + O. Sauter, C. Angioni and Y. R. Lin-Liu, + Physics of Plasmas 6 (1999) 2834 + O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) + Physics of Plasmas 9 (2002) 5140 + """ + zz = zef[j - 1] + zft = tpf(j, triang, sqeps) + _nues = nues(j, rmajor, zef, mu, sqeps, tempe, ne) + _nues_sqrt = np.sqrt(_nues) + zdf = (1.0 + 0.26 * (1.0 - zft) * _nues_sqrt) + ( + 0.18 * (1.0 - 0.37 * zft) * _nues / np.sqrt(zz) + ) + zfte = zft / zdf # $f^{32\_ee}_{teff}(\nu_{e*})$, Eq.15d + zfte2 = zfte * zfte + zfte3 = zfte * zfte2 + zfte4 = zfte2 * zfte2 + + zdf = (1.0 + (1.0 + 0.6 * zft) * _nues_sqrt) + ( + 0.85 * (1.0 - 0.37 * zft) * _nues * (1.0 + zz) + ) + zfti = zft / zdf # $f^{32\_ei}_{teff}(\nu_{e*})$, Eq.15e + zfti2 = zfti * zfti + zfti3 = zfti * zfti2 + zfti4 = zfti2 * zfti2 + + # $F_{32\_ee}(X)$, Eq.15b + hcee = ( + ((0.05 + 0.62 * zz) / zz / (1.0 + 0.44 * zz) * (zfte - zfte4)) + + ((zfte2 - zfte4 - 1.2 * (zfte3 - zfte4)) / (1.0 + 0.22 * zz)) + + (1.2 / (1.0 + 0.5 * zz) * zfte4) + ) + + # $F_{32\_ei}(Y)$, Eq.15c + hcei = ( + (-(0.56 + 1.93 * zz) / zz / (1.0 + 0.44 * zz) * (zfti - zfti4)) + + (4.95 / (1.0 + 2.48 * zz) * (zfti2 - zfti4 - 0.55 * (zfti3 - zfti4))) + - (1.2 / (1.0 + 0.5 * zz) * zfti4) + ) + + # Corrections suggested by Fable, 15/05/2015 + return beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) * ( + hcee + hcei + ) + dcsa( + j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps + ) * beta_poloidal_local( + j, nr, rmajor, bt, ne, tempe, mu, rho + ) / beta_poloidal_local_total( + j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho + ) + + +@nb.jit(nopython=True, cache=True) +def xcsa( + j, + nr, + rmajor, + bt, + triang, + mu, + sqeps, + tempi, + tempe, + amain, + zmain, + ni, + ne, + rho, + zef, +): + """Grad(ln(Ti)) coefficient in the Sauter bootstrap scaling + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 2 to nr + nr : input integer : maximum value of j + This function calculates the coefficient scaling grad(ln(Ti)) + in the Sauter bootstrap current scaling. + Code by Angioni, 29th May 2002. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + O. Sauter, C. Angioni and Y. R. Lin-Liu, + Physics of Plasmas 6 (1999) 2834 + O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) + Physics of Plasmas 9 (2002) 5140 + """ + zz = zef[j - 1] + zft = tpf(j, triang, sqeps) + _nues = nues(j, rmajor, zef, mu, sqeps, tempe, ne) + zdf = (1.0 + (1.0 - 0.1 * zft) * np.sqrt(_nues)) + 0.5 * ( + 1.0 - 0.5 * zft + ) * _nues / zz + zfte = zft / zdf # $f^{34}_{teff}(\nu_{e*})$, Eq.16b + + # Eq.16a + xcsa = ((1.0 + 1.4 / (zz + 1.0)) * zfte - 1.9 / (zz + 1.0) * zfte * zfte) + ( + 0.3 * zfte * zfte + 0.2 * zfte * zfte * zfte + ) * zfte / (zz + 1.0) + + # $\alpha_0$, Eq.17a + a0 = (-1.17 * (1.0 - zft)) / (1.0 - 0.22 * zft - 0.19 * zft * zft) + + _nuis = nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni) + _nuis_sqrt = np.sqrt(_nuis) + a1 = _nuis**2 * zft**6 + + # $\alpha(\nu_{i*})$, Eq.17b + alp = ( + (a0 + 0.25 * (1.0 - zft * zft) * _nuis_sqrt) / (1.0 + 0.5 * _nuis_sqrt) + + 0.315 * a1 + ) / (1.0 + 0.15 * a1) + + # Corrections suggested by Fable, 15/05/2015 + return ( + beta_poloidal_local_total(j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho) + - beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) + ) * (xcsa * alp) + dcsa( + j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps + ) * ( + 1.0 + - beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) + / beta_poloidal_local_total(j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho) + ) + + +@nb.jit(nopython=True, cache=True) +def nues(j, rmajor, zef, mu, sqeps, tempe, ne): + """Relative frequency of electron collisions + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + This function calculates the relative frequency of electron + collisions: NU* = Nuei*q*Rt/eps**1.5/Vte + The electron-ion collision frequency NUEI=NUEE*1.4*ZEF is + used. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + Yushmanov, 30th April 1987 (?) + """ + return ( + nuee(j, tempe, ne) + * 1.4 + * zef[j - 1] + * rmajor + / np.abs(mu[j - 1] * (sqeps[j - 1] ** 3) * np.sqrt(tempe[j - 1]) * 1.875e7) + ) + + +@nb.jit(nopython=True, cache=True) +def beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho): + """Local beta poloidal calculation + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + nr : input integer : maximum value of j + This function calculates the local beta poloidal. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). +

beta poloidal = 4*pi*ne*Te/Bpo**2 + Pereverzev, 25th April 1989 (?) + """ + return ( + np.where( + j != nr, + 1.6e-4 * np.pi * (ne[j] + ne[j - 1]) * (tempe[j] + tempe[j - 1]), + 6.4e-4 * np.pi * ne[j - 1] * tempe[j - 1], + ) + * (rmajor / (bt * rho[j - 1] * np.abs(mu[j - 1] + 1.0e-4))) ** 2 + ) + + +@nb.jit(nopython=True, cache=True) +def beta_poloidal_local_total(j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho): + """Local beta poloidal calculation, including ion pressure + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + nr : input integer : maximum value of j + This function calculates the local total beta poloidal. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). +

beta poloidal = 4*pi*(ne*Te+ni*Ti)/Bpo**2 + where ni is the sum of all ion densities (thermal) + Pereverzev, 25th April 1989 (?) + E Fable, private communication, 15th May 2014 + """ + + return ( + np.where( + j != nr, + 1.6e-4 + * np.pi + * ( + ((ne[j] + ne[j - 1]) * (tempe[j] + tempe[j - 1])) + + ((ni[j] + ni[j - 1]) * (tempi[j] + tempi[j - 1])) + ), + 6.4e-4 * np.pi * (ne[j - 1] * tempe[j - 1] + ni[j - 1] * tempi[j - 1]), + ) + * (rmajor / (bt * rho[j - 1] * np.abs(mu[j - 1] + 1.0e-4))) ** 2 + ) + + +@nb.jit(nopython=True, cache=True) +def tpf(j, triang, sqeps, fit=1): + """Trapped particle fraction + author: P J Knight, CCFE, Culham Science Centre + j : input integer : radial element index in range 1 to nr + fit : input integer : (1)=ASTRA method, 2=Equation from Sauter2002, 3=Equation from Sauter2013 + This function calculates the trapped particle fraction at + a given radius. +

A number of different fits are provided, but the one + to be used is hardwired prior to run-time. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + O. Sauter et al, Plasma Phys. Contr. Fusion 44 (2002) 1999 + O. Sauter, 2013: + http://infoscience.epfl.ch/record/187521/files/lrp_012013.pdf + """ + s = sqeps[j - 1] + eps = s * s + + if fit == 1: + # ASTRA method, from Emiliano Fable, private communication + # (Excluding h term which dominates for inverse aspect ratios < 0.5, + # and tends to take the trapped particle fraction to 1) + + zz = 1.0 - eps + return 1.0 - zz * np.sqrt(zz) / (1.0 + 1.46 * s) + elif fit == 2: + # Equation 4 of Sauter 2002 + # Similar to, but not quite identical to g above + + return 1.0 - (1.0 - eps) ** 2 / (1.0 + 1.46 * s) / np.sqrt(1.0 - eps * eps) + elif fit == 3: + # Includes correction for triangularity + + epseff = 0.67 * (1.0 - 1.4 * triang * np.abs(triang)) * eps + + return 1.0 - np.sqrt((1.0 - eps) / (1.0 + eps)) * (1.0 - epseff) / ( + 1.0 + 2.0 * np.sqrt(epseff) + ) + + raise RuntimeError(f"fit={fit} is not valid. Must be 1, 2, or 3") + + class Physics: def __init__(self, plasma_profile, current_drive): self.outfile = constants.nout @@ -234,19 +982,17 @@ def physics(self): ) # Hender scaling for diamagnetic current at tight physics_variables.aspect ratio - current_drive_variables.diacf_hender = self.diamagnetic_fraction_hender( + current_drive_variables.diacf_hender = diamagnetic_fraction_hender( physics_variables.beta ) # SCENE scaling for diamagnetic current - current_drive_variables.diacf_scene = self.diamagnetic_fraction_scene( + current_drive_variables.diacf_scene = diamagnetic_fraction_scene( physics_variables.beta, physics_variables.q95, physics_variables.q0 ) # Pfirsch-Schlüter scaling for diamagnetic current - current_drive_variables.pscf_scene = self.ps_fraction_scene( - physics_variables.beta - ) + current_drive_variables.pscf_scene = ps_fraction_scene(physics_variables.beta) current_drive_variables.bscf_sauter = ( current_drive_variables.cboot * self.bootstrap_fraction_sauter() @@ -426,7 +1172,7 @@ def physics(self): # Calculate ion/electron equilibration power - physics_variables.piepv = self.rether( + physics_variables.piepv = rether( physics_variables.alphan, physics_variables.alphat, physics_variables.dene, @@ -599,7 +1345,7 @@ def physics(self): physics_variables.vsind, physics_variables.vsres, physics_variables.vsstt, - ) = self.vscalc( + ) = vscalc( physics_variables.csawth, physics_variables.eps, physics_variables.facoh, @@ -611,6 +1357,7 @@ def physics(self): times_variables.t_fusion_ramp, times_variables.tburn, physics_variables.rli, + constants.rmu0, ) # Calculate auxiliary physics related information @@ -666,7 +1413,7 @@ def physics(self): # Hartmann and Zohm physics_variables.dnbeta = 4.0e0 * physics_variables.rli - physics_variables.betalim = self.culblm( + physics_variables.betalim = culblm( physics_variables.bt, physics_variables.dnbeta, physics_variables.plascur, @@ -821,12 +1568,15 @@ def physics(self): po.write( self.outfile, ( - f" 'fzactual, frac, reinke_variables.impvardiv = {reinke_variables.fzactual}, {impurity_radiation_module.impurity_arr_frac(reinke_variables.impvardiv)}, {reinke_variables.impvardiv}" + f" 'fzactual, frac, reinke_variables.impvardiv = {reinke_variables.fzactual}," + f" {impurity_radiation_module.impurity_arr_frac(reinke_variables.impvardiv)}," + f" {reinke_variables.impvardiv}" ), ) + @staticmethod def culdlm( - self, bt, idensl, pdivt, plascur, prn1, qcyl, q95, rmajor, rminor, sarea, zeff + bt, idensl, pdivt, plascur, prn1, qcyl, q95, rmajor, rminor, sarea, zeff ): """Density limit calculation author: P J Knight, CCFE, Culham Science Centre @@ -865,16 +1615,18 @@ def culdlm( # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. - dlim = 1.54e20 * qperp**0.43 * bt**0.31 / (q95 * rmajor) ** 0.45 - dlimit[0] = dlim / prn1 + dlimit[0] = ( + 1.54e20 * qperp**0.43 * bt**0.31 / (q95 * rmajor) ** 0.45 + ) / prn1 # Borrass density limit model for ITER (I) # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. # Borrass et al, ITER-TN-PH-9-6 (1989) - dlim = 1.8e20 * qperp**0.53 * bt**0.31 / (q95 * rmajor) ** 0.22 - dlimit[1] = dlim / prn1 + dlimit[1] = ( + 1.8e20 * qperp**0.53 * bt**0.31 / (q95 * rmajor) ** 0.22 + ) / prn1 # Borrass density limit model for ITER (II) # This applies to the density at the plasma edge, so must be scaled @@ -882,8 +1634,9 @@ def culdlm( # This formula is (almost) identical to that in the original routine # denlim (now deleted). - dlim = 0.5e20 * qperp**0.57 * bt**0.31 / (q95 * rmajor) ** 0.09 - dlimit[2] = dlim / prn1 + dlimit[2] = ( + 0.5e20 * qperp**0.57 * bt**0.31 / (q95 * rmajor) ** 0.09 + ) / prn1 # JET edge radiation density limit model # This applies to the density at the plasma edge, so must be scaled @@ -900,15 +1653,13 @@ def culdlm( dlimit[3] = 0.0 else: - dlim = 1.0e20 * np.sqrt(pdivt / denom) - dlimit[3] = dlim / prn1 + dlimit[3] = (1.0e20 * np.sqrt(pdivt / denom)) / prn1 # JET simplified density limit model # This applies to the density at the plasma edge, so must be scaled # to give the density limit applying to the average plasma density. - dlim = 0.237e20 * bt * np.sqrt(pdivt) / rmajor - dlimit[4] = dlim / prn1 + dlimit[4] = (0.237e20 * bt * np.sqrt(pdivt) / rmajor) / prn1 # Hugill-Murakami M.q limit # qcyl=qstar here, which is okay according to the literature @@ -923,7 +1674,8 @@ def culdlm( return dlimit, dlimit[idensl - 1] - def plasma_composition(self): + @staticmethod + def plasma_composition(): """Calculates various plasma component fractional makeups author: P J Knight, CCFE, Culham Science Centre @@ -1131,8 +1883,8 @@ def plasma_composition(self): / impurity_radiation_module.impurity_arr_amass[imp] ) + @staticmethod def phyaux( - self, aspect, dene, deni, @@ -1210,30 +1962,8 @@ def phyaux( return burnup, dntau, figmer, fusrat, qfuel, rndfuel, taup - def rether(self, alphan, alphat, dene, dlamie, te, ti, zeffai): - """Routine to find the equilibration power between the - ions and electrons - author: P J Knight, CCFE, Culham Science Centre - alphan : input real : density profile index - alphat : input real : temperature profile index - dene : input real : electron density (/m3) - dlamie : input real : ion-electron coulomb logarithm - te : input real : electron temperature (keV) - ti : input real : ion temperature (keV) - zeffai : input real : mass weighted plasma effective charge - piepv : output real : ion/electron equilibration power (MW/m3) - This routine calculates the equilibration power between the - ions and electrons. - Unknown origin - """ - profie = (1.0 + alphan) ** 2 / ( - (2.0 * alphan - 0.5 * alphat + 1.0) * np.sqrt(1.0 + alphat) - ) - conie = 2.42165e-41 * dlamie * dene**2 * zeffai * profie - - return conie * (ti - te) / (te**1.5) - - def pohm(self, facoh, kappa95, plascur, rmajor, rminor, ten, vol, zeff): + @staticmethod + def pohm(facoh, kappa95, plascur, rmajor, rminor, ten, vol, zeff): # Density weighted electron temperature in 10 keV units t10 = ten / 10.0 @@ -1274,117 +2004,8 @@ def pohm(self, facoh, kappa95, plascur, rmajor, rminor, ten, vol, zeff): return pohmpv, pohmmw, rpfac, rplas - def vscalc( - self, - csawth, - eps, - facoh, - gamma, - kappa, - rmajor, - rplas, - plascur, - t_fusion_ramp, - tburn, - rli, - ): - """Volt-second requirements - author: P J Knight, CCFE, Culham Science Centre - csawth : input real : coefficient for sawteeth effects - eps : input real : inverse aspect ratio - facoh : input real : fraction of plasma current produced inductively - gamma : input real : Ejima coeff for resistive start-up V-s component - kappa : input real : plasma elongation - plascur: input real : plasma current (A) - rli : input real : plasma normalised inductivity - rmajor : input real : plasma major radius (m) - rplas : input real : plasma resistance (ohm) - t_fusion_ramp : input real : heating time (s) - tburn : input real : burn time (s) - phiint : output real : internal plasma volt-seconds (Wb) - rlp : output real : plasma inductance (H) - vsbrn : output real : volt-seconds needed during flat-top (heat+burn) (Wb) - vsind : output real : internal and external plasma inductance V-s (Wb) - vsres : output real : resistive losses in start-up volt-seconds (Wb) - vsstt : output real : total volt-seconds needed (Wb) - This subroutine calculates the volt-second requirements and some - other related items. - AEA FUS 251: A User's Guide to the PROCESS Systems Code - """ - # Internal inductance - - rlpint = constants.rmu0 * rmajor * rli / 2.0 - phiint = rlpint * plascur - - # Start-up resistive component - # Uses ITER formula without the 10 V-s add-on - - vsres = gamma * constants.rmu0 * plascur * rmajor - - # Hirshman, Neilson: Physics of Fluids, 29 (1986) p790 - # fit for external inductance - - aeps = (1.0 + 1.81 * np.sqrt(eps) + 2.05 * eps) * np.log(8.0 / eps) - ( - 2.0 + 9.25 * np.sqrt(eps) - 1.21 * eps - ) - beps = ( - 0.73 - * np.sqrt(eps) - * (1.0 + 2.0 * eps**4 - 6.0 * eps**5 + 3.7 * eps**6) - ) - rlpext = ( - rmajor * constants.rmu0 * aeps * (1.0 - eps) / (1.0 - eps + beps * kappa) - ) - - rlp = rlpext + rlpint - - # Inductive V-s component - - vsind = rlp * plascur - vsstt = vsres + vsind - - # Loop voltage during flat-top - # Include enhancement factor in flattop V-s requirement - # to account for MHD sawtooth effects. - - vburn = plascur * rplas * facoh * csawth - - # N.B. tburn on first iteration will not be correct - # if the pulsed reactor option is used, but the value - # will be correct on subsequent calls. - - vsbrn = vburn * (t_fusion_ramp + tburn) - vsstt = vsstt + vsbrn - - return phiint, rlp, vsbrn, vsind, vsres, vsstt - - def culblm(self, bt, dnbeta, plascur, rminor): - """Beta scaling limit - author: P J Knight, CCFE, Culham Science Centre - bt : input real : toroidal B-field on plasma axis (T) - dnbeta : input real : Troyon-like g coefficient - plascur : input real : plasma current (A) - rminor : input real : plasma minor axis (m) - betalim : output real : beta limit as defined below - This subroutine calculates the beta limit, using - the algorithm documented in AEA FUS 172. - The limit applies to beta defined with respect to the total B-field. - Switch iculbl determines which components of beta to include. - - If iculbl = 0, then the limit is applied to the total beta - If iculbl = 1, then the limit is applied to the thermal beta only - If iculbl = 2, then the limit is applied to the thermal + neutral beam beta components - If iculbl = 3, then the limit is applied to the toroidal beta - - The default value for the g coefficient is dnbeta = 3.5 - AEA FUS 172: Physics Assessment for the European Reactor Study - AEA FUS 251: A User's Guide to the PROCESS Systems Code - """ - - return 0.01 * dnbeta * (plascur / 1.0e6) / (rminor * bt) - + @staticmethod def culcur( - self, alphaj, alphap, bt, @@ -1473,9 +2094,7 @@ def culcur( if icurr == 1: # Peng analytical fit fq = (1.22 - 0.68 * eps) / ((1.0 - eps * eps) ** 2) * sf**2 elif icurr == 2: # Peng scaling for double null divertor; TARTs [STAR Code] - curhat = ( - 1.0e6 * physics_module.plasc(qpsi, asp, rminor, bt, kappa, triang) / bt - ) + curhat = 1.0e6 * plasc(qpsi, asp, rminor, bt, kappa, triang) / bt elif icurr == 3: # Simple ITER scaling (simply the cylindrical case) fq = 1.0 elif icurr == 4: # ITER formula (IPDG89) @@ -1504,7 +2123,7 @@ def culcur( fq *= 1 if icurr == 7 else (1.0 + (abs(kappa95 - 1.2)) ** 3) elif icurr == 7: # Connor-Hastie asymptotically-correct expression # N.B. If iprofile=1, alphaj will be wrong during the first call (only) - fq = physics_module.conhas(alphaj, alphap, bt, triang95, eps, kappa95, p0) + fq = conhas(alphaj, alphap, bt, triang95, eps, kappa95, p0, constants.rmu0) elif ( icurr == 8 ): # Sauter scaling allowing negative triangularity [FED May 2016] @@ -1542,7 +2161,7 @@ def culcur( physics_variables.normalised_total_beta = ( 1.0e8 * physics_variables.beta * rminor * bt / plascur ) - bp = self.bpol(icurr, plascur, qpsi, asp, bt, kappa, triang, pperim) + bp = bpol(icurr, plascur, qpsi, asp, 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... @@ -1553,2545 +2172,1550 @@ def culcur( return alphaj, rli, bp, qstar, plascur - def conhas(self, alphaj, alphap, bt, delta95, eps, kappa95, p0): - """Routine to calculate the F coefficient used for scaling the - plasma current - author: P J Knight, CCFE, Culham Science Centre - alphaj : input real : current profile index - alphap : input real : pressure profile index - bt : input real : toroidal field on axis (T) - delta95 : input real : plasma triangularity 95% - eps : input real : inverse aspect ratio - kappa95 : input real : plasma elongation 95% - p0 : input real : central plasma pressure (Pa) - fq : output real : scaling for edge q from circular - cross-section cylindrical case - This routine calculates the F coefficient used for scaling the - plasma current, using the Connor-Hastie scaling given in - AEA FUS 172. - AEA FUS 172: Physics Assessment for the European Reactor Study - AEA FUS 251: A User's Guide to the PROCESS Systems Code - """ - - # Exponent in Connor-Hastie current profile - matching total - # current gives the following trivial relation - - lamda = alphaj - - # Exponent in Connor-Hastie pressure profile - - nu = alphap - - # Central plasma beta - - beta0 = 2.0 * constants.rmu0 * p0 / (bt**2) - - # Plasma internal inductance - - lamp1 = 1.0 + lamda - li = lamp1 / lamda * (lamp1 / lamda * np.log(lamp1) - 1.0) - - # T/r in AEA FUS 172 - - kap1 = kappa95 + 1.0 - tr = kappa95 * delta95 / kap1**2 - - # E/r in AEA FUS 172 - - er = (kappa95 - 1.0) / kap1 + @staticmethod + def eped_warning(): + eped_warning = "" - # T primed in AEA FUS 172 + if (physics_variables.triang < 0.399e0) or (physics_variables.triang > 0.601e0): + eped_warning += f"{physics_variables.triang = }" - tprime = 2.0 * tr * lamp1 / (1.0 + 0.5 * lamda) + if (physics_variables.kappa < 1.499e0) or (physics_variables.kappa > 2.001e0): + eped_warning += f"{physics_variables.kappa = }" - # E primed in AEA FUS 172 + if (physics_variables.plascur < 9.99e6) or ( + physics_variables.plascur > 20.01e6 + ): + eped_warning += f"{physics_variables.plascur = }" - eprime = er * lamp1 / (1.0 + lamda / 3.0) + if (physics_variables.rmajor < 6.99e0) or (physics_variables.rmajor > 11.01e0): + eped_warning += f"{physics_variables.rmajor = }" - # Delta primed in AEA FUS 172 + if (physics_variables.rminor < 1.99e0) or (physics_variables.rminor > 3.501e0): + eped_warning += f"{physics_variables.rminor = }" - deltap = 0.5 * kap1 * eps * 0.5 * li + beta0 / ( - 0.5 * kap1 * eps - ) * lamp1**2 / (1.0 + nu) + if (physics_variables.normalised_total_beta < 1.99e0) or ( + physics_variables.normalised_total_beta > 3.01e0 + ): + eped_warning += f"{physics_variables.normalised_total_beta = }" - # Delta/R0 in AEA FUS 172 + if physics_variables.tesep > 0.5: + eped_warning += f"{physics_variables.tesep = }" - deltar = beta0 / 6.0 * (1.0 + 5.0 * lamda / 6.0 + 0.25 * lamda**2) + ( - 0.5 * kap1 * eps - ) ** 2 * 0.125 * (1.0 - (lamda**2) / 3.0) + return eped_warning - # F coefficient + def outtim(self): + po.oheadr(self.outfile, "Times") - return (0.5 * kap1) ** 2 * ( - 1.0 - + eps**2 * (0.5 * kap1) ** 2 - + 0.5 * deltap**2 - + 2.0 * deltar - + 0.5 * (eprime**2 + er**2) - + 0.5 * (tprime**2 + 4.0 * tr**2) + po.ovarrf( + self.outfile, + "Initial charge time for CS from zero current (s)", + "(tramp)", + times_variables.tramp, + ) + po.ovarrf( + self.outfile, + "Plasma current ramp-up time (s)", + "(tohs)", + times_variables.tohs, + ) + po.ovarrf( + self.outfile, + "Heating time (s)", + "(t_fusion_ramp)", + times_variables.t_fusion_ramp, + ) + po.ovarre( + self.outfile, "Burn time (s)", "(tburn)", times_variables.tburn, "OP " + ) + po.ovarrf( + self.outfile, + "Reset time to zero current for CS (s)", + "(tqnch)", + times_variables.tqnch, + ) + po.ovarrf( + self.outfile, "Time between pulses (s)", "(tdwell)", times_variables.tdwell + ) + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Total plant cycle time (s)", + "(tcycle)", + times_variables.tcycle, + "OP ", ) - def plasc(self, qbar, aspect, rminor, bt, kappa, delta): - """Function to calculate plasma current (Peng scaling) - author: J Galambos, FEDC/ORNL + def outplas(self): + """Subroutine to output the plasma physics information author: P J Knight, CCFE, Culham Science Centre - aspect : input real : plasma aspect ratio - bt : input real : toroidal field on axis (T) - delta : input real : plasma triangularity - kappa : input real : plasma elongation - qbar : input real : edge q-bar - rminor : input real : plasma minor radius (m) - This function calculates the plasma current in MA, - using a scaling from Peng, Galambos and Shipe (1992). - It is primarily used for Tight Aspect Ratio Tokamaks and is - selected via icurr=2. - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, - Fusion Technology, 21, 1729 + self.outfile : input integer : Fortran output unit identifier + This routine writes the plasma physics information + to a file, in a tidy format. + AEA FUS 251: A User's Guide to the PROCESS Systems Code """ - eps = 1.0 / aspect - - c1 = kappa**2 / (1.0 + delta) + delta - c2 = kappa**2 / (1.0 - delta) - delta - - d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 - d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 - - if aspect < c1: - y1 = np.sqrt((c1 * eps - 1.0) / (1.0 + eps)) * (1.0 + delta) / kappa - else: - y1 = np.sqrt((1.0 - c1 * eps) / (1.0 + eps)) * (1.0 + delta) / kappa - - y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * (1.0 - delta) / kappa - - e1 = 2.0 * kappa / (d1 * (1.0 + delta)) - e2 = 2.0 * kappa / (d2 * (1.0 - delta)) - - h2 = (1.0 + (c2 - 1.0) * eps / 2.0) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) - f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * (c2 * eps + 1.0)) - g = eps * kappa / (1.0 - eps * delta) - ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) - - if aspect < c1: - h1 = (1.0 + (1.0 - c1) * eps / 2.0) / np.sqrt( - (1.0 + eps) * (c1 * eps - 1.0) + # ############################################### + # Dimensionless plasma parameters. See reference below. + physics_module.nu_star = ( + 1 + / constants.rmu0 + * (15.0e0 * constants.echarge**4 * physics_variables.dlamie) + / (4.0e0 * np.pi**1.5e0 * constants.epsilon0**2) + * physics_variables.vol**2 + * physics_variables.rmajor**2 + * physics_variables.bt + * np.sqrt(physics_variables.eps) + * physics_variables.dnla**3 + * physics_variables.kappa + / ( + physics_module.total_plasma_internal_energy**2 + * physics_variables.plascur ) - f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) - ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) - else: - h1 = (1.0 + (1.0 - c1) * eps / 2.0) / np.sqrt( - (1.0 + eps) * (1.0 - c1 * eps) - ) - f1 = -(d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) - ff1 = f1 * (-g + 2.0 * h1 * np.arctan(y1)) - - return ( - rminor - * bt - / qbar - * 5.0 - * kappa - / (2.0 * np.pi**2) - * (np.arcsin(e1) / e1 + np.arcsin(e2) / e2) - * (ff1 + ff2) - ) - - def bpol(self, icurr, ip, qbar, aspect, bt, kappa, delta, perim): - """Function to calculate poloidal field - author: J Galambos, FEDC/ORNL - author: P J Knight, CCFE, Culham Science Centre - icurr : input integer : current scaling model to use - ip : input real : plasma current (A) - qbar : input real : edge q-bar - aspect : input real : plasma aspect ratio - bt : input real : toroidal field on axis (T) - kappa : input real : plasma elongation - delta : input real : plasma triangularity - perim : input real : plasma perimeter (m) - This function calculates the poloidal field in Tesla, - using a simple calculation using Stoke's Law for conventional - tokamaks, or for TARTs, a scaling from Peng, Galambos and - Shipe (1992). - J D Galambos, STAR Code : Spherical Tokamak Analysis and Reactor Code, - unpublished internal Oak Ridge document - Y.-K. M. Peng, J. Galambos and P.C. Shipe, 1992, - Fusion Technology, 21, 1729 - """ - if icurr != 2: - return constants.rmu0 * ip / perim - - # Original coding, only suitable for TARTs [STAR Code] + ) - eps = 1.0 / aspect + physics_module.rho_star = np.sqrt( + 2.0e0 + * constants.mproton + * physics_variables.aion + * physics_module.total_plasma_internal_energy + / (3.0e0 * physics_variables.vol * physics_variables.dnla) + ) / ( + constants.echarge + * physics_variables.bt + * physics_variables.eps + * physics_variables.rmajor + ) - c1 = kappa**2 / (1.0 + delta) + delta - c2 = kappa**2 / (1.0 - delta) - delta + physics_module.beta_mcdonald = ( + 4.0e0 + / 3.0e0 + * constants.rmu0 + * physics_module.total_plasma_internal_energy + / (physics_variables.vol * physics_variables.bt**2) + ) - d1 = (kappa / (1.0 + delta)) ** 2 + 1.0 - d2 = (kappa / (1.0 - delta)) ** 2 + 1.0 + po.oheadr(self.outfile, "Plasma") - if aspect < c1: - y1 = np.sqrt((c1 * eps - 1.0) / (1.0 + eps)) * (1.0 + delta) / kappa + if stellarator_variables.istell == 0: + if physics_variables.idivrt == 0: + po.ocmmnt(self.outfile, "Plasma configuration = limiter") + elif physics_variables.idivrt == 1: + po.ocmmnt(self.outfile, "Plasma configuration = single null divertor") + elif physics_variables.idivrt == 2: + po.ocmmnt(self.outfile, "Plasma configuration = double null divertor") + else: + error_handling.idiags[0] = physics_variables.idivrt + po.report_error(85) else: - y1 = np.sqrt((1.0 - c1 * eps) / (1.0 + eps)) * (1.0 + delta) / kappa - - y2 = np.sqrt((c2 * eps + 1.0) / (1.0 - eps)) * (1.0 - delta) / kappa - - h2 = (1.0 + (c2 - 1.0) * eps / 2.0) / np.sqrt((1.0 - eps) * (c2 * eps + 1.0)) - f2 = (d2 * (1.0 - delta) * eps) / ((1.0 - eps) * (c2 * eps + 1.0)) - g = eps * kappa / (1.0 - eps * delta) - ff2 = f2 * (g + 2.0 * h2 * np.arctan(y2)) + po.ocmmnt(self.outfile, "Plasma configuration = stellarator") - if aspect < c1: - h1 = (1.0 + (1.0 - c1) * eps / 2.0) / np.sqrt( - (1.0 + eps) * (c1 * eps - 1.0) - ) - f1 = (d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) - ff1 = f1 * (g - h1 * np.log((1.0 + y1) / (1.0 - y1))) - else: - h1 = (1.0 + (1.0 - c1) * eps / 2.0) / np.sqrt( - (1.0 + eps) * (1.0 - c1 * eps) - ) - f1 = -(d1 * (1.0 + delta) * eps) / ((1.0 + eps) * (c1 * eps - 1.0)) - ff1 = f1 * (-g + 2.0 * h1 * np.arctan(y1)) + if stellarator_variables.istell == 0: + if physics_variables.itart == 0: + physics_module.itart_r = physics_variables.itart + po.ovarrf( + self.outfile, + "Tokamak aspect ratio = Conventional, itart = 0", + "(itart)", + physics_module.itart_r, + ) + elif physics_variables.itart == 1: + physics_module.itart_r = physics_variables.itart + po.ovarrf( + self.outfile, + "Tokamak aspect ratio = Spherical, itart = 1", + "(itart)", + physics_module.itart_r, + ) - return bt * (ff1 + ff2) / (2.0 * np.pi * qbar) + po.osubhd(self.outfile, "Plasma Geometry :") + po.ovarrf( + self.outfile, "Major radius (m)", "(rmajor)", physics_variables.rmajor + ) + po.ovarrf( + self.outfile, + "Minor radius (m)", + "(rminor)", + physics_variables.rminor, + "OP ", + ) + po.ovarrf(self.outfile, "Aspect ratio", "(aspect)", physics_variables.aspect) - def bootstrap_fraction_iter89( - self, aspect, beta, bt, cboot, plascur, q95, q0, rmajor, vol - ): - """Original ITER calculation of bootstrap-driven fraction - of the plasma current. - author: P J Knight, CCFE, Culham Science Centre - aspect : input real : plasma aspect ratio - beta : input real : plasma total beta - bt : input real : toroidal field on axis (T) - cboot : input real : bootstrap current fraction multiplier - plascur : input real : plasma current (A) - q95 : input real : safety factor at 95% surface - q0 : input real : central safety factor - rmajor : input real : plasma major radius (m) - vol : input real : plasma volume (m3) - This routine performs the original ITER calculation of the - plasma current bootstrap fraction. - ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, - ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 - """ - xbs = min(10, q95 / q0) - cbs = cboot * (1.32 - 0.235 * xbs + 0.0185 * xbs**2) - bpbs = ( - constants.rmu0 - * plascur - / (2 * np.pi * np.sqrt(vol / (2 * np.pi**2 * rmajor))) - ) - betapbs = beta * bt**2 / bpbs**2 + if stellarator_variables.istell == 0: + if physics_variables.ishape in [0, 6, 8]: + po.ovarrf( + self.outfile, + "Elongation, X-point (input value used)", + "(kappa)", + physics_variables.kappa, + "IP ", + ) + elif physics_variables.ishape == 1: + po.ovarrf( + self.outfile, + "Elongation, X-point (TART scaling)", + "(kappa)", + physics_variables.kappa, + "OP ", + ) + elif physics_variables.ishape in [2, 3]: + po.ovarrf( + self.outfile, + "Elongation, X-point (Zohm scaling)", + "(kappa)", + physics_variables.kappa, + "OP ", + ) + po.ovarrf( + self.outfile, + "Zohm scaling adjustment factor", + "(fkzohm)", + physics_variables.fkzohm, + ) + elif physics_variables.ishape in [4, 5, 7]: + po.ovarrf( + self.outfile, + "Elongation, X-point (calculated from kappa95)", + "(kappa)", + physics_variables.kappa, + "OP ", + ) + elif physics_variables.ishape == 9: + po.ovarrf( + self.outfile, + "Elongation, X-point (calculated from aspect ratio and li(3))", + "(kappa)", + physics_variables.kappa, + "OP ", + ) + elif physics_variables.ishape == 10: + po.ovarrf( + self.outfile, + "Elongation, X-point (calculated from aspect ratio and stability margin)", + "(kappa)", + physics_variables.kappa, + "OP ", + ) + elif physics_variables.ishape == 11: + po.ovarrf( + self.outfile, + "Elongation, X-point (calculated from aspect ratio via Menard 2016)", + "(kappa)", + physics_variables.kappa, + "OP ", + ) + else: + error_handling.idiags[0] = physics_variables.ishape + po.report_error(86) - if betapbs <= 0.0: # only possible if beta <= 0.0 - return 0.0 - return cbs * (betapbs / np.sqrt(aspect)) ** 1.3 + if physics_variables.ishape in [4, 5, 7]: + po.ovarrf( + self.outfile, + "Elongation, 95% surface (input value used)", + "(kappa95)", + physics_variables.kappa95, + "IP ", + ) + else: + po.ovarrf( + self.outfile, + "Elongation, 95% surface (calculated from kappa)", + "(kappa95)", + physics_variables.kappa95, + "OP ", + ) - def bootstrap_fraction_nevins( - self, - alphan, - alphat, - betat, - bt, - dene, - plascur, - q95, - q0, - rmajor, - rminor, - ten, - zeff, - ): - """Bootstrap current fraction from Nevins et al scaling - author: P J Knight, CCFE, Culham Science Centre - alphan : input real : density profile index - alphat : input real : temperature profile index - betat : input real : total plasma beta (with respect to the toroidal - field) - bt : input real : toroidal field on axis (T) - dene : input real : electron density (/m3) - plascur: input real : plasma current (A) - q0 : input real : central safety factor - q95 : input real : safety factor at 95% surface - rmajor : input real : plasma major radius (m) - rminor : input real : plasma minor radius (m) - ten : input real : density weighted average plasma temperature (keV) - zeff : input real : plasma effective charge - This function calculates the bootstrap current fraction, - using the Nevins et al method, 4/11/90. - AEA FUS 251: A User's Guide to the PROCESS Systems Code - """ - # Calculate peak electron beta + po.ovarrf( + self.outfile, + "Elongation, area ratio calc.", + "(kappaa)", + physics_variables.kappaa, + "OP ", + ) - betae0 = ( - physics_variables.ne0 - * physics_variables.te0 - * 1.0e3 - * constants.echarge - / (bt**2 / (2.0 * constants.rmu0)) - ) + if physics_variables.ishape in [0, 2, 6, 8, 9, 10, 11]: + po.ovarrf( + self.outfile, + "Triangularity, X-point (input value used)", + "(triang)", + physics_variables.triang, + "IP ", + ) + elif physics_variables.ishape == 1: + po.ovarrf( + self.outfile, + "Triangularity, X-point (TART scaling)", + "(triang)", + physics_variables.triang, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Triangularity, X-point (calculated from triang95)", + "(triang)", + physics_variables.triang, + "OP ", + ) - # Call integration routine + if physics_variables.ishape in [3, 4, 5, 7]: + po.ovarrf( + self.outfile, + "Triangularity, 95% surface (input value used)", + "(triang95)", + physics_variables.triang95, + "IP ", + ) + else: + po.ovarrf( + self.outfile, + "Triangularity, 95% surface (calculated from triang)", + "(triang95)", + physics_variables.triang95, + "OP ", + ) - ainteg, _ = integrate.quad( - lambda y: self.bsinteg( - y, dene, ten, bt, rminor, rmajor, zeff, alphat, alphan, q0, q95, betat - ), - 0, - 0.999, - epsabs=0.001, - epsrel=0.001, - ) + po.ovarrf( + self.outfile, + "Plasma poloidal perimeter (m)", + "(pperim)", + physics_variables.pperim, + "OP ", + ) - # Calculate bootstrap current and fraction + po.ovarrf( + self.outfile, + "Plasma cross-sectional area (m2)", + "(xarea)", + physics_variables.xarea, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma surface area (m2)", + "(sarea)", + physics_variables.sarea, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma volume (m3)", + "(vol)", + physics_variables.vol, + "OP ", + ) - aibs = 2.5 * betae0 * rmajor * bt * q95 * ainteg - return 1.0e6 * aibs / plascur + po.osubhd(self.outfile, "Current and Field :") - def bsinteg( - self, y, dene, ten, bt, rminor, rmajor, zeff, alphat, alphan, q0, q95, betat - ): - """Integrand function for Nevins et al bootstrap current scaling - author: P J Knight, CCFE, Culham Science Centre - y : input real : abscissa of integration, = normalised minor radius - This function calculates the integrand function for the - Nevins et al bootstrap current scaling, 4/11/90. - AEA FUS 251: A User's Guide to the PROCESS Systems Code - """ - # Constants for fit to q-profile + if stellarator_variables.istell == 0: + if physics_variables.iprofile == 0: + po.ocmmnt( + self.outfile, + "Consistency between q0,q,alphaj,rli,dnbeta is not enforced", + ) + else: + po.ocmmnt( + self.outfile, + "Consistency between q0,q,alphaj,rli,dnbeta is enforced", + ) - c1 = 1.0 - c2 = 1.0 - c3 = 1.0 + po.oblnkl(self.outfile) + po.ovarin( + self.outfile, + "Plasma current scaling law used", + "(icurr)", + physics_variables.icurr, + ) - # Compute average electron beta - - betae = ( - dene * ten * 1.0e3 * constants.echarge / (bt**2 / (2.0 * constants.rmu0)) - ) - - nabla = rminor * np.sqrt(y) / rmajor - x = (1.46 * np.sqrt(nabla) + 2.4 * nabla) / (1.0 - nabla) ** 1.5 - z = zeff - d = ( - 1.414 * z - + z * z - + x * (0.754 + 2.657 * z + 2.0 * z * z) - + x * x * (0.348 + 1.243 * z + z * z) - ) - al2 = -x * (0.884 + 2.074 * z) / d - a2 = alphat * (1.0 - y) ** (alphan + alphat - 1.0) - alphai = -1.172 / (1.0 + 0.462 * x) - a1 = (alphan + alphat) * (1.0 - y) ** (alphan + alphat - 1.0) - al1 = x * (0.754 + 2.21 * z + z * z + x * (0.348 + 1.243 * z + z * z)) / d - - # q-profile - - q = q0 + (q95 - q0) * (c1 * y + c2 * y * y + c3 * y**3) / (c1 + c2 + c3) - - pratio = (betat - betae) / betae - - return (q / q95) * (al1 * (a1 + pratio * (a1 + alphai * a2)) + al2 * a2) - - def bootstrap_fraction_wilson( - self, alphaj, alphap, alphat, betpth, q0, qpsi, rmajor, rminor - ): - """Bootstrap current fraction from Wilson et al scaling - author: P J Knight, CCFE, Culham Science Centre - alphaj : input real : current profile index - alphap : input real : pressure profile index - alphat : input real : temperature profile index - beta : input real : total beta - betpth : input real : thermal component of poloidal beta - q0 : input real : safety factor on axis - qpsi : input real : edge safety factor - rmajor : input real : major radius (m) - rminor : input real : minor radius (m) - This function calculates the bootstrap current fraction - using the numerically fitted algorithm written by Howard Wilson. - AEA FUS 172: Physics Assessment for the European Reactor Study - H. R. Wilson, Nuclear Fusion 32 (1992) 257 - """ - term1 = np.log(0.5) - term2 = np.log(q0 / qpsi) + po.ovarrf( + self.outfile, + "Plasma current (MA)", + "(plascur/1D6)", + physics_variables.plascur / 1.0e6, + "OP ", + ) - termp = 1.0 - 0.5 ** (1.0 / alphap) - termt = 1.0 - 0.5 ** (1.0 / alphat) - termj = 1.0 - 0.5 ** (1.0 / alphaj) + if physics_variables.iprofile == 1: + po.ovarrf( + self.outfile, + "Current density profile factor", + "(alphaj)", + physics_variables.alphaj, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Current density profile factor", + "(alphaj)", + physics_variables.alphaj, + ) - alfpnw = term1 / np.log(np.log((q0 + (qpsi - q0) * termp) / qpsi) / term2) - alftnw = term1 / np.log(np.log((q0 + (qpsi - q0) * termt) / qpsi) / term2) - aj = term1 / np.log(np.log((q0 + (qpsi - q0) * termj) / qpsi) / term2) + po.ovarrf( + self.outfile, + "Plasma internal inductance, li", + "(rli)", + physics_variables.rli, + "OP ", + ) + po.ovarrf( + self.outfile, + "Vertical field at plasma (T)", + "(bvert)", + physics_variables.bvert, + "OP ", + ) - # Crude check for NaN errors or other illegal values... + po.ovarrf( + self.outfile, + "Vacuum toroidal field at R (T)", + "(bt)", + physics_variables.bt, + ) + po.ovarrf( + self.outfile, + "Average poloidal field (T)", + "(bp)", + physics_variables.bp, + "OP ", + ) - if np.isnan(aj) or np.isnan(alfpnw) or np.isnan(alftnw) or aj < 0: - error_handling.fdiags[0] = aj - error_handling.fdiags[1] = alfpnw - error_handling.fdiags[2] = alftnw - error_handling.fdiags[3] = aj + po.ovarrf( + self.outfile, + "Total field (sqrt(bp^2 + bt^2)) (T)", + "(btot)", + physics_variables.btot, + "OP ", + ) - error_handling.report_error(76) + if stellarator_variables.istell == 0: + po.ovarrf( + self.outfile, "Safety factor on axis", "(q0)", physics_variables.q0 + ) - # Ratio of ionic charge to electron charge + if physics_variables.icurr == 2: + po.ovarrf( + self.outfile, "Mean edge safety factor", "(q)", physics_variables.q + ) - z = 1.0 + po.ovarrf( + self.outfile, + "Safety factor at 95% flux surface", + "(q95)", + physics_variables.q95, + ) - # Inverse aspect ratio: r2 = maximum plasma radius, r1 = minimum + po.ovarrf( + self.outfile, + "Cylindrical safety factor (qcyl)", + "(qstar)", + physics_variables.qstar, + "OP ", + ) - r2 = rmajor + rminor - r1 = rmajor - rminor - eps1 = (r2 - r1) / (r2 + r1) + if physics_variables.ishape == 1: + po.ovarrf( + self.outfile, + "Lower limit for edge safety factor q", + "(qlim)", + physics_variables.qlim, + "OP ", + ) - # Coefficients fitted using least squares techniques + else: + po.ovarrf( + self.outfile, + "Rotational transform", + "(iotabar)", + stellarator_variables.iotabar, + ) - saj = np.sqrt(aj) + po.osubhd(self.outfile, "Beta Information :") - a = np.array( - [ - 1.41 * (1.0 - 0.28 * saj) * (1.0 + 0.12 / z), - 0.36 * (1.0 - 0.59 * saj) * (1.0 + 0.8 / z), - -0.27 * (1.0 - 0.47 * saj) * (1.0 + 3.0 / z), - 0.0053 * (1.0 + 5.0 / z), - -0.93 * (1.0 - 0.34 * saj) * (1.0 + 0.15 / z), - -0.26 * (1.0 - 0.57 * saj) * (1.0 - 0.27 * z), - 0.064 * (1.0 - 0.6 * aj + 0.15 * aj * aj) * (1.0 + 7.6 / z), - -0.0011 * (1.0 + 9.0 / z), - -0.33 * (1.0 - aj + 0.33 * aj * aj), - -0.26 * (1.0 - 0.87 / saj - 0.16 * aj), - -0.14 * (1.0 - 1.14 / saj - 0.45 * saj), - -0.0069, - ] + betath = ( + physics_variables.beta - physics_variables.betaft - physics_variables.betanb ) + gammaft = (physics_variables.betaft + physics_variables.betanb) / betath - seps1 = np.sqrt(eps1) - - b = np.array( - [ - 1.0, - alfpnw, - alftnw, - alfpnw * alftnw, - seps1, - alfpnw * seps1, - alftnw * seps1, - alfpnw * alftnw * seps1, - eps1, - alfpnw * eps1, - alftnw * eps1, - alfpnw * alftnw * eps1, - ] + po.ovarre(self.outfile, "Total plasma beta", "(beta)", physics_variables.beta) + po.ovarre( + self.outfile, + "Total poloidal beta", + "(betap)", + physics_variables.betap, + "OP ", + ) + po.ovarre( + self.outfile, + "Total toroidal beta", + " ", + physics_variables.beta + * (physics_variables.btot / physics_variables.bt) ** 2, + "OP ", + ) + po.ovarre( + self.outfile, "Fast alpha beta", "(betaft)", physics_variables.betaft, "OP " + ) + po.ovarre( + self.outfile, "Beam ion beta", "(betanb)", physics_variables.betanb, "OP " + ) + po.ovarre( + self.outfile, + "(Fast alpha + beam physics_variables.beta)/(thermal physics_variables.beta)", + "(gammaft)", + gammaft, + "OP ", ) - # Empirical bootstrap current fraction - - return seps1 * betpth * (a * b).sum() + po.ovarre(self.outfile, "Thermal beta", " ", betath, "OP ") + po.ovarre( + self.outfile, + "Thermal poloidal beta", + " ", + betath * (physics_variables.btot / physics_variables.bp) ** 2, + "OP ", + ) + po.ovarre( + self.outfile, + "Thermal toroidal physics_variables.beta (= beta-exp)", + " ", + betath * (physics_variables.btot / physics_variables.bt) ** 2, + "OP ", + ) - def diamagnetic_fraction_hender(self, beta): - """author: S.I. Muldrew, CCFE, Culham Science Centre - Diamagnetic contribution at tight aspect ratio. - Tim Hender fit - """ - return beta / 2.8 + po.ovarrf( + self.outfile, + "2nd stability physics_variables.beta : beta_p / (R/a)", + "(eps*betap)", + physics_variables.eps * physics_variables.betap, + "OP ", + ) + po.ovarrf( + self.outfile, + "2nd stability physics_variables.beta upper limit", + "(epbetmax)", + physics_variables.epbetmax, + ) - def diamagnetic_fraction_scene(self, beta, q95, q0): - """author: S.I. Muldrew, CCFE, Culham Science Centre - Diamagnetic fraction based on SCENE fit by Tim Hender - See Issue #992 - """ - return beta * (0.1 * q95 / q0 + 0.44) * 4.14e-1 + if stellarator_variables.istell == 0: + if physics_variables.iprofile == 1: + po.ovarrf( + self.outfile, + "Beta g coefficient", + "(dnbeta)", + physics_variables.dnbeta, + "OP ", + ) + else: + po.ovarrf( + self.outfile, + "Beta g coefficient", + "(dnbeta)", + physics_variables.dnbeta, + ) - def ps_fraction_scene(self, beta): - """author: S.I. Muldrew, CCFE, Culham Science Centre - Pfirsch-Schlüter fraction based on SCENE fit by Tim Hender - See Issue #992 - """ - return -9e-2 * beta - - def bootstrap_fraction_sauter(self): - """Bootstrap current fraction from Sauter et al scaling - author: P J Knight, CCFE, Culham Science Centre - None - This function calculates the bootstrap current fraction - using the Sauter, Angioni and Lin-Liu scaling. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 - """ - NR = 200 - - roa = np.arange(1, NR + 1, step=1) / NR + po.ovarrf( + self.outfile, + "Normalised thermal beta", + " ", + 1.0e8 + * betath + * physics_variables.rminor + * physics_variables.bt + / physics_variables.plascur, + "OP ", + ) - rho = np.sqrt(physics_variables.xarea / np.pi) * roa - sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) + po.ovarrf( + self.outfile, + "Normalised total beta", + " ", + physics_variables.normalised_total_beta, + "OP ", + ) - ne = 1e-19 * np.vectorize( - lambda r: profiles_module.nprofile( - r, - physics_variables.rhopedn, - physics_variables.ne0, - physics_variables.neped, - physics_variables.nesep, - physics_variables.alphan, + normalised_toroidal_beta = ( + physics_variables.normalised_total_beta + * (physics_variables.btot / physics_variables.bt) ** 2 ) - )(roa) - ni = (physics_variables.dnitot / physics_variables.dene) * ne - tempe = np.vectorize( - lambda r: profiles_module.tprofile( - r, - physics_variables.rhopedt, - physics_variables.te0, - physics_variables.teped, - physics_variables.tesep, - physics_variables.alphat, - physics_variables.tbeta, + po.ovarrf( + self.outfile, + "Normalised toroidal beta", + "(normalised_toroidal_beta)", + normalised_toroidal_beta, + "OP ", ) - )(roa) - tempi = (physics_variables.ti / physics_variables.te) * tempe - - zef = np.full_like(tempi, physics_variables.zeff) # Flat Zeff profile assumed - - # mu = 1/safety factor - # Parabolic q profile assumed - mu = 1 / ( - physics_variables.q0 - + (physics_variables.q - physics_variables.q0) * roa**2 - ) - amain = np.full_like(mu, physics_variables.afuel) - zmain = np.full_like(mu, 1.0 + physics_variables.fhe3) - - if ne[NR - 1] == 0.0: - ne[NR - 1] = 1e-4 * ne[NR - 2] - ni[NR - 1] = 1e-4 * ni[NR - 2] - - if tempe[NR - 1] == 0.0: - tempe[NR - 1] = 1e-4 * tempe[NR - 2] - tempi[NR - 1] = 1e-4 * tempi[NR - 2] - - # Calculate total bootstrap current (MA) by summing along profiles - iboot = 0.0 - for ir in range(0, NR): - if ir + 1 == NR: - jboot = 0.0 - da = 0.0 - else: - drho = rho[ir + 1] - rho[ir] - da = 2 * np.pi * rho[ir] * drho # area of annulus - - dlogte_drho = (np.log(tempe[ir + 1]) - np.log(tempe[ir])) / drho - dlogti_drho = (np.log(tempi[ir + 1]) - np.log(tempi[ir])) / drho - dlogne_drho = (np.log(ne[ir + 1]) - np.log(ne[ir])) / drho - - # The factor of 0.5 below arises because in ASTRA the logarithms - # are coded as (e.g.): (Te(j+1)-Te(j))/(Te(j+1)+Te(j)), which - # actually corresponds to grad(log(Te))/2. So the factors dcsa etc. - # are a factor two larger than one might otherwise expect. - jboot = 0.5 * ( - self.dcsa( - ir + 1, - NR, - physics_variables.rmajor, - physics_variables.bt, - physics_variables.triang, - ne, - ni, - tempe, - tempi, - mu, - rho, - zef, - sqeps, - ) - * dlogne_drho - + self.hcsa( - ir + 1, - NR, - physics_variables.rmajor, - physics_variables.bt, - physics_variables.triang, - ne, - ni, - tempe, - tempi, - mu, - rho, - zef, - sqeps, - ) - * dlogte_drho - + self.xcsa( - ir + 1, - NR, - physics_variables.rmajor, - physics_variables.bt, - physics_variables.triang, - mu, - sqeps, - tempi, - tempe, - amain, - zmain, - ni, - ne, - rho, - zef, - ) - * dlogti_drho - ) - jboot = ( - -physics_variables.bt - / (0.2 * np.pi * physics_variables.rmajor) - * rho[ir] - * mu[ir] - * jboot - ) # MA/m2 - - iboot += da * jboot - return 1.0e6 * iboot / physics_variables.plascur - - def hcsa( - self, j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps - ): - """Grad(ln(Te)) coefficient in the Sauter bootstrap scaling - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - nr : input integer : maximum value of j - This function calculates the coefficient scaling grad(ln(Te)) - in the Sauter bootstrap current scaling. - Code by Angioni, 29th May 2002. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 - """ - - if j == 1: - return 0.0 - - zz = zef[j - 1] - zft = self.tpf(j, triang, sqeps) - zdf = 1.0 + 0.26 * (1.0 - zft) * np.sqrt( - self.nues(j, rmajor, zef, mu, sqeps, tempe, ne) - ) - zdf = zdf + 0.18 * (1.0 - 0.37 * zft) * self.nues( - j, rmajor, zef, mu, sqeps, tempe, ne - ) / np.sqrt(zz) - zfte = zft / zdf # $f^{32\_ee}_{teff}(\nu_{e*})$, Eq.15d - zfte2 = zfte * zfte - zfte3 = zfte * zfte2 - zfte4 = zfte2 * zfte2 - - zdf = 1.0 + (1.0 + 0.6 * zft) * np.sqrt( - self.nues(j, rmajor, zef, mu, sqeps, tempe, ne) - ) - zdf = zdf + 0.85 * (1.0 - 0.37 * zft) * self.nues( - j, rmajor, zef, mu, sqeps, tempe, ne - ) * (1.0 + zz) - zfti = zft / zdf # $f^{32\_ei}_{teff}(\nu_{e*})$, Eq.15e - zfti2 = zfti * zfti - zfti3 = zfti * zfti2 - zfti4 = zfti2 * zfti2 - - hcee = (0.05 + 0.62 * zz) / zz / (1.0 + 0.44 * zz) * (zfte - zfte4) - hcee = hcee + (zfte2 - zfte4 - 1.2 * (zfte3 - zfte4)) / (1.0 + 0.22 * zz) - hcee = hcee + 1.2 / (1.0 + 0.5 * zz) * zfte4 # $F_{32\_ee}(X)$, Eq.15b - - hcei = -(0.56 + 1.93 * zz) / zz / (1.0 + 0.44 * zz) * (zfti - zfti4) - hcei = hcei + 4.95 / (1.0 + 2.48 * zz) * ( - zfti2 - zfti4 - 0.55 * (zfti3 - zfti4) - ) - hcei = hcei - 1.2 / (1.0 + 0.5 * zz) * zfti4 # $F_{32\_ei}(Y)$, Eq.15c - - # Corrections suggested by Fable, 15/05/2015 - return self.beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) * ( - hcee + hcei - ) + self.dcsa( - j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps - ) * self.beta_poloidal_local( - j, nr, rmajor, bt, ne, tempe, mu, rho - ) / self.beta_poloidal_local_total( - j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho - ) - - def xcsa( - self, - j, - nr, - rmajor, - bt, - triang, - mu, - sqeps, - tempi, - tempe, - amain, - zmain, - ni, - ne, - rho, - zef, - ): - """Grad(ln(Ti)) coefficient in the Sauter bootstrap scaling - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - nr : input integer : maximum value of j - This function calculates the coefficient scaling grad(ln(Ti)) - in the Sauter bootstrap current scaling. - Code by Angioni, 29th May 2002. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 - """ - if j == 1: - return 0.0 - - zz = zef[j - 1] - zft = self.tpf(j, triang, sqeps) - zdf = 1.0 + (1.0 - 0.1 * zft) * np.sqrt( - self.nues(j, rmajor, zef, mu, sqeps, tempe, ne) - ) - zdf = ( - zdf - + 0.5 - * (1.0 - 0.5 * zft) - * self.nues(j, rmajor, zef, mu, sqeps, tempe, ne) - / zz - ) - zfte = zft / zdf # $f^{34}_{teff}(\nu_{e*})$, Eq.16b - - xcsa = (1.0 + 1.4 / (zz + 1.0)) * zfte - 1.9 / (zz + 1.0) * zfte * zfte - xcsa = xcsa + (0.3 * zfte * zfte + 0.2 * zfte * zfte * zfte) * zfte / ( - zz + 1.0 - ) # Eq.16a - - a0 = -1.17 * (1.0 - zft) - a0 = a0 / (1.0 - 0.22 * zft - 0.19 * zft * zft) # $\alpha_0$, Eq.17a - - alp = ( - a0 - + 0.25 - * (1.0 - zft * zft) - * np.sqrt(self.nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni)) - ) / ( - 1.0 - + 0.5 * np.sqrt(self.nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni)) - ) - a1 = ( - self.nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni) - * self.nuis(j, rmajor, mu, sqeps, tempi, amain, zmain, ni) - * zft**6 - ) - alp = (alp + 0.315 * a1) / (1.0 + 0.15 * a1) # $\alpha(\nu_{i*})$, Eq.17b - - # Corrections suggested by Fable, 15/05/2015 - return ( - self.beta_poloidal_local_total( - j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho + if physics_variables.iculbl == 0: + po.ovarrf( + self.outfile, + "Limit on total beta", + "(betalim)", + physics_variables.betalim, + "OP ", ) - - self.beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) - ) * (xcsa * alp) + self.dcsa( - j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps - ) * ( - 1.0 - - self.beta_poloidal_local(j, nr, rmajor, bt, ne, tempe, mu, rho) - / self.beta_poloidal_local_total( - j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho + elif physics_variables.iculbl == 1: + po.ovarrf( + self.outfile, + "Limit on thermal beta", + "(betalim)", + physics_variables.betalim, + "OP ", ) - ) - - def nuis(self, j, rmajor, mu, sqeps, tempi, amain, zmain, ni): - """Relative frequency of ion collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the relative frequency of ion - collisions: NU* = Nui*q*Rt/eps**1.5/Vti - The full ion collision frequency NUI is used. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - Yushmanov, 30th April 1987 (?) - """ - return ( - 3.2e-6 - * self.nui(j, zmain, ni, tempi, amain) - * rmajor - / ( - abs(mu[j - 1] + 1.0e-4) - * sqeps[j - 1] ** 3 - * np.sqrt(tempi[j - 1] / amain[j - 1]) + else: + po.ovarrf( + self.outfile, + "Limit on thermal + NB beta", + "(betalim)", + physics_variables.betalim, + "OP ", ) - ) - def nui(self, j, zmain, ni, tempi, amain): - """Full frequency of ion collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the full frequency of ion - collisions (Hz). -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - None - """ - # Coulomb logarithm = 15 is used - return ( - zmain[j - 1] ** 4 - * ni[j - 1] - * 322.0 - / (tempi[j - 1] * np.sqrt(tempi[j - 1] * amain[j - 1])) + po.ovarre( + self.outfile, + "Plasma thermal energy (J)", + " ", + 1.5e0 + * betath + * physics_variables.btot + * physics_variables.btot + / (2.0e0 * constants.rmu0) + * physics_variables.vol, + "OP ", ) - def dcsa( - self, j, nr, rmajor, bt, triang, ne, ni, tempe, tempi, mu, rho, zef, sqeps - ): - """Grad(ln(ne)) coefficient in the Sauter bootstrap scaling - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - nr : input integer : maximum value of j - This function calculates the coefficient scaling grad(ln(ne)) - in the Sauter bootstrap current scaling. - Code by Angioni, 29th May 2002. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter, C. Angioni and Y. R. Lin-Liu, - Physics of Plasmas 6 (1999) 2834 - O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) - Physics of Plasmas 9 (2002) 5140 - - DCSA $\\equiv \\mathcal{L}_{31}$, Eq.14a, Sauter et al, 1999 - """ - - if j == 1: - return 0.0 + po.ovarre( + self.outfile, + "Total plasma internal energy (J)", + "(total_plasma_internal_energy)", + physics_module.total_plasma_internal_energy, + "OP ", + ) - zz = zef[j - 1] - zft = self.tpf(j, triang, sqeps) - zdf = 1.0 + (1.0 - 0.1 * zft) * np.sqrt( - self.nues(j, rmajor, zef, mu, sqeps, tempe, ne) + po.osubhd(self.outfile, "Temperature and Density (volume averaged) :") + po.ovarrf( + self.outfile, "Electron temperature (keV)", "(te)", physics_variables.te + ) + po.ovarrf( + self.outfile, + "Electron temperature on axis (keV)", + "(te0)", + physics_variables.te0, + "OP ", ) - zdf = ( - zdf - + 0.5 * (1.0 - zft) * self.nues(j, rmajor, zef, mu, sqeps, tempe, ne) / zz + po.ovarrf(self.outfile, "Ion temperature (keV)", "(ti)", physics_variables.ti) + po.ovarrf( + self.outfile, + "Ion temperature on axis (keV)", + "(ti0)", + physics_variables.ti0, + "OP ", ) - zft = zft / zdf # $f^{31}_{teff}(\nu_{e*})$, Eq.14b - dcsa = ( - (1.0 + 1.4 / (zz + 1.0)) * zft - - 1.9 / (zz + 1.0) * zft * zft - + (0.3 * zft * zft + 0.2 * zft * zft * zft) * zft / (zz + 1.0) + po.ovarrf( + self.outfile, + "Electron temp., density weighted (keV)", + "(ten)", + physics_variables.ten, + "OP ", ) - - # Corrections suggested by Fable, 15/05/2015 - return dcsa * self.beta_poloidal_local_total( - j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho + po.ovarre( + self.outfile, "Electron density (/m3)", "(dene)", physics_variables.dene ) - - def nues(self, j, rmajor, zef, mu, sqeps, tempe, ne): - """Relative frequency of electron collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the relative frequency of electron - collisions: NU* = Nuei*q*Rt/eps**1.5/Vte - The electron-ion collision frequency NUEI=NUEE*1.4*ZEF is - used. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - Yushmanov, 30th April 1987 (?) - """ - return ( - self.nuee(j, tempe, ne) - * 1.4 - * zef[j - 1] - * rmajor - / abs(mu[j - 1] * (sqeps[j - 1] ** 3) * np.sqrt(tempe[j - 1]) * 1.875e7) + po.ovarre( + self.outfile, + "Electron density on axis (/m3)", + "(ne0)", + physics_variables.ne0, + "OP ", ) - - def nuee(self, j, tempe, ne): - """Frequency of electron-electron collisions - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the frequency of electron-electron - collisions (Hz): NUEE = 4*SQRT(pi)/3*Ne*e**4*lambd/ - SQRT(Me)/Te**1.5 -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - Yushmanov, 25th April 1987 (?), - updated by Pereverzev, 9th November 1994 (?) - """ - return ( - 670.0 - * self.coulg(j, tempe, ne) - * ne[j - 1] - / (tempe[j - 1] * np.sqrt(tempe[j - 1])) + po.ovarre( + self.outfile, + "Line-averaged electron density (/m3)", + "(dnla)", + physics_variables.dnla, + "OP ", ) - def coulg(self, j, tempe, ne): - """Coulomb logarithm - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - This function calculates the Coulomb logarithm, valid - for e-e collisions (T_e > 0.01 keV), and for - e-i collisions (T_e > 0.01*Zeff^2) (Alexander, 9/5/1994). -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - C. A. Ordonez and M. I. Molina, Phys. Plasmas 1 (1994) 2515 - Rev. Mod. Phys., V.48, Part 1 (1976) 275 - """ - return 15.9 - 0.5 * np.log(ne[j - 1]) + np.log(tempe[j - 1]) - - def beta_poloidal_local(self, j, nr, rmajor, bt, ne, tempe, mu, rho): - """Local beta poloidal calculation - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - nr : input integer : maximum value of j - This function calculates the local beta poloidal. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). -

beta poloidal = 4*pi*ne*Te/Bpo**2 - Pereverzev, 25th April 1989 (?) - """ - if j != nr: - beta_poloidal_local = ( - 1.6e-4 * np.pi * (ne[j] + ne[j - 1]) * (tempe[j] + tempe[j - 1]) + if stellarator_variables.istell == 0: + po.ovarre( + self.outfile, + "Line-averaged electron density / Greenwald density", + "(dnla_gw)", + physics_variables.dnla / physics_variables.dlimit[6], + "OP ", ) - else: - beta_poloidal_local = 6.4e-4 * np.pi * ne[j - 1] * tempe[j - 1] - - return ( - beta_poloidal_local - * (rmajor / (bt * rho[j - 1] * abs(mu[j - 1] + 1.0e-4))) ** 2 - ) - def beta_poloidal_local_total( - self, j, nr, rmajor, bt, ne, ni, tempe, tempi, mu, rho - ): - """Local beta poloidal calculation, including ion pressure - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - nr : input integer : maximum value of j - This function calculates the local total beta poloidal. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). -

beta poloidal = 4*pi*(ne*Te+ni*Ti)/Bpo**2 - where ni is the sum of all ion densities (thermal) - Pereverzev, 25th April 1989 (?) - E Fable, private communication, 15th May 2014 - """ - - if j != nr: - beta_poloidal_local_total = ( - 1.6e-4 - * np.pi - * ( - ((ne[j] + ne[j - 1]) * (tempe[j] + tempe[j - 1])) - + ((ni[j] + ni[j - 1]) * (tempi[j] + tempi[j - 1])) - ) - ) - else: - beta_poloidal_local_total = ( - 6.4e-4 * np.pi * (ne[j - 1] * tempe[j - 1] + ni[j - 1] * tempi[j - 1]) - ) - - return ( - beta_poloidal_local_total - * (rmajor / (bt * rho[j - 1] * abs(mu[j - 1] + 1.0e-4))) ** 2 - ) - - def tpf(self, j, triang, sqeps, fit=1): - """Trapped particle fraction - author: P J Knight, CCFE, Culham Science Centre - j : input integer : radial element index in range 1 to nr - fit : input integer : (1)=ASTRA method, 2=Equation from Sauter2002, 3=Equation from Sauter2013 - This function calculates the trapped particle fraction at - a given radius. -

A number of different fits are provided, but the one - to be used is hardwired prior to run-time. -

The code was supplied by Emiliano Fable, IPP Garching - (private communication). - O. Sauter et al, Plasma Phys. Contr. Fusion 44 (2002) 1999 - O. Sauter, 2013: - http://infoscience.epfl.ch/record/187521/files/lrp_012013.pdf - """ - s = sqeps[j - 1] - eps = s * s - - if fit == 1: - # ASTRA method, from Emiliano Fable, private communication - # (Excluding h term which dominates for inverse aspect ratios < 0.5, - # and tends to take the trapped particle fraction to 1) - - zz = 1.0 - eps - return 1.0 - zz * np.sqrt(zz) / (1.0 + 1.46 * s) - elif fit == 2: - # Equation 4 of Sauter 2002 - # Similar to, but not quite identical to g above - - return 1.0 - (1.0 - eps) ** 2 / (1.0 + 1.46 * s) / np.sqrt(1.0 - eps * eps) - elif fit == 3: - # Includes correction for triangularity - - epseff = 0.67 * (1.0 - 1.4 * triang * abs(triang)) * eps - - return 1.0 - np.sqrt((1.0 - eps) / (1.0 + eps)) * (1.0 - epseff) / ( - 1.0 + 2.0 * np.sqrt(epseff) - ) - - raise RuntimeError(f"{fit=} is not valid. Must be 1, 2, or 3") - - def eped_warning(self): - eped_warning = "" - - if (physics_variables.triang < 0.399e0) or (physics_variables.triang > 0.601e0): - eped_warning += f"{physics_variables.triang = }" - - if (physics_variables.kappa < 1.499e0) or (physics_variables.kappa > 2.001e0): - eped_warning += f"{physics_variables.kappa = }" - - if (physics_variables.plascur < 9.99e6) or ( - physics_variables.plascur > 20.01e6 - ): - eped_warning += f"{physics_variables.plascur = }" - - if (physics_variables.rmajor < 6.99e0) or (physics_variables.rmajor > 11.01e0): - eped_warning += f"{physics_variables.rmajor = }" - - if (physics_variables.rminor < 1.99e0) or (physics_variables.rminor > 3.501e0): - eped_warning += f"{physics_variables.rminor = }" - - if (physics_variables.normalised_total_beta < 1.99e0) or ( - physics_variables.normalised_total_beta > 3.01e0 - ): - eped_warning += f"{physics_variables.normalised_total_beta = }" - - if physics_variables.tesep > 0.5: - eped_warning += f"{physics_variables.tesep = }" - - return eped_warning - - def outtim(self): - po.oheadr(self.outfile, "Times") - - po.ovarrf( + po.ovarre( self.outfile, - "Initial charge time for CS from zero current (s)", - "(tramp)", - times_variables.tramp, + "Ion density (/m3)", + "(dnitot)", + physics_variables.dnitot, + "OP ", ) - po.ovarrf( - self.outfile, - "Plasma current ramp-up time (s)", - "(tohs)", - times_variables.tohs, + po.ovarre( + self.outfile, "Fuel density (/m3)", "(deni)", physics_variables.deni, "OP " ) - po.ovarrf( + po.ovarre( self.outfile, - "Heating time (s)", - "(t_fusion_ramp)", - times_variables.t_fusion_ramp, + "Total impurity density with Z > 2 (no He) (/m3)", + "(dnz)", + physics_variables.dnz, + "OP ", ) po.ovarre( - self.outfile, "Burn time (s)", "(tburn)", times_variables.tburn, "OP " - ) - po.ovarrf( self.outfile, - "Reset time to zero current for CS (s)", - "(tqnch)", - times_variables.tqnch, - ) - po.ovarrf( - self.outfile, "Time between pulses (s)", "(tdwell)", times_variables.tdwell + "Helium ion density (thermalised ions only) (/m3)", + "(dnalp)", + physics_variables.dnalp, + "OP ", ) - po.oblnkl(self.outfile) po.ovarre( self.outfile, - "Total plant cycle time (s)", - "(tcycle)", - times_variables.tcycle, + "Proton density (/m3)", + "(dnprot)", + physics_variables.dnprot, "OP ", ) - - def outplas(self): - """Subroutine to output the plasma physics information - author: P J Knight, CCFE, Culham Science Centre - self.outfile : input integer : Fortran output unit identifier - This routine writes the plasma physics information - to a file, in a tidy format. - AEA FUS 251: A User's Guide to the PROCESS Systems Code - """ - - # ############################################### - # Dimensionless plasma parameters. See reference below. - physics_module.nu_star = ( - 1 - / constants.rmu0 - * (15.0e0 * constants.echarge**4 * physics_variables.dlamie) - / (4.0e0 * np.pi**1.5e0 * constants.epsilon0**2) - * physics_variables.vol**2 - * physics_variables.rmajor**2 - * physics_variables.bt - * np.sqrt(physics_variables.eps) - * physics_variables.dnla**3 - * physics_variables.kappa - / ( - physics_module.total_plasma_internal_energy**2 - * physics_variables.plascur + if physics_variables.protium > 1.0e-10: + po.ovarre( + self.outfile, + "Seeded protium density / electron density", + "(protium)", + physics_variables.protium, ) - ) - physics_module.rho_star = np.sqrt( - 2.0e0 - * constants.mproton - * physics_variables.aion - * physics_module.total_plasma_internal_energy - / (3.0e0 * physics_variables.vol * physics_variables.dnla) - ) / ( - constants.echarge - * physics_variables.bt - * physics_variables.eps - * physics_variables.rmajor + po.ovarre( + self.outfile, + "Hot beam density (/m3)", + "(dnbeam)", + physics_variables.dnbeam, + "OP ", ) - - physics_module.beta_mcdonald = ( - 4.0e0 - / 3.0e0 - * constants.rmu0 - * physics_module.total_plasma_internal_energy - / (physics_variables.vol * physics_variables.bt**2) + po.ovarre( + self.outfile, + "Density limit from scaling (/m3)", + "(dnelimt)", + physics_variables.dnelimt, + "OP ", ) + if (numerics.ioptimz > 0) and (numerics.active_constraints[4]): + po.ovarre( + self.outfile, + "Density limit (enforced) (/m3)", + "(boundu(9)*dnelimt)", + numerics.boundu[8] * physics_variables.dnelimt, + "OP ", + ) - po.oheadr(self.outfile, "Plasma") + po.ovarre( + self.outfile, + "Helium ion density (thermalised ions only) / electron density", + "(ralpne)", + physics_variables.ralpne, + ) + po.oblnkl(self.outfile) - if stellarator_variables.istell == 0: - if physics_variables.idivrt == 0: - po.ocmmnt(self.outfile, "Plasma configuration = limiter") - elif physics_variables.idivrt == 1: - po.ocmmnt(self.outfile, "Plasma configuration = single null divertor") - elif physics_variables.idivrt == 2: - po.ocmmnt(self.outfile, "Plasma configuration = double null divertor") - else: - error_handling.idiags[0] = physics_variables.idivrt - po.report_error(85) - else: - po.ocmmnt(self.outfile, "Plasma configuration = stellarator") + po.ocmmnt(self.outfile, "Impurities") + po.oblnkl(self.outfile) + po.ocmmnt(self.outfile, "Plasma ion densities / electron density:") - if stellarator_variables.istell == 0: - if physics_variables.itart == 0: - physics_module.itart_r = physics_variables.itart - po.ovarrf( - self.outfile, - "Tokamak aspect ratio = Conventional, itart = 0", - "(itart)", - physics_module.itart_r, + for imp in range(impurity_radiation_module.nimp): + # MDK Update fimp, as this will make the ITV output work correctly. + impurity_radiation_module.fimp[ + imp + ] = impurity_radiation_module.impurity_arr_frac[imp] + str1 = ( + f2py_compatible_to_string( + impurity_radiation_module.impurity_arr_label[imp] ) - elif physics_variables.itart == 1: - physics_module.itart_r = physics_variables.itart - po.ovarrf( - self.outfile, - "Tokamak aspect ratio = Spherical, itart = 1", - "(itart)", - physics_module.itart_r, + + " concentration" + ) + str2 = f"(fimp({imp + 1:02}))" + # MDK Add output flag for H which is calculated. + if imp == 0: + po.ovarre( + self.outfile, str1, str2, impurity_radiation_module.fimp[imp], "OP " ) + else: + po.ovarre(self.outfile, str1, str2, impurity_radiation_module.fimp[imp]) - po.osubhd(self.outfile, "Plasma Geometry :") - po.ovarrf( - self.outfile, "Major radius (m)", "(rmajor)", physics_variables.rmajor - ) - po.ovarrf( + po.ovarre( self.outfile, - "Minor radius (m)", - "(rminor)", - physics_variables.rminor, + "Average mass of all ions (amu)", + "(aion)", + physics_variables.aion, "OP ", ) - po.ovarrf(self.outfile, "Aspect ratio", "(aspect)", physics_variables.aspect) + # MDK Say which impurity is varied, if iteration variable fimpvar (102) is turned on + # if (any(ixc == 102)) : + # call ovarst(self.outfile,'Impurity used as an iteration variable' , '', '"' // impurity_arr(impvar)%label // '"') + # po.ovarre(self.outfile,'Fractional density of variable impurity (ion / electron density)','(fimpvar)',fimpvar) + # + po.oblnkl(self.outfile) + po.ovarrf( + self.outfile, "Effective charge", "(zeff)", physics_variables.zeff, "OP " + ) - if stellarator_variables.istell == 0: - if physics_variables.ishape in [0, 6, 8]: - po.ovarrf( - self.outfile, - "Elongation, X-point (input value used)", - "(kappa)", - physics_variables.kappa, - "IP ", - ) - elif physics_variables.ishape == 1: - po.ovarrf( - self.outfile, - "Elongation, X-point (TART scaling)", - "(kappa)", - physics_variables.kappa, - "OP ", - ) - elif physics_variables.ishape in [2, 3]: - po.ovarrf( - self.outfile, - "Elongation, X-point (Zohm scaling)", - "(kappa)", - physics_variables.kappa, - "OP ", - ) - po.ovarrf( - self.outfile, - "Zohm scaling adjustment factor", - "(fkzohm)", - physics_variables.fkzohm, - ) - elif physics_variables.ishape in [4, 5, 7]: - po.ovarrf( - self.outfile, - "Elongation, X-point (calculated from kappa95)", - "(kappa)", - physics_variables.kappa, - "OP ", - ) - elif physics_variables.ishape == 9: - po.ovarrf( - self.outfile, - "Elongation, X-point (calculated from aspect ratio and li(3))", - "(kappa)", - physics_variables.kappa, - "OP ", - ) - elif physics_variables.ishape == 10: - po.ovarrf( - self.outfile, - "Elongation, X-point (calculated from aspect ratio and stability margin)", - "(kappa)", - physics_variables.kappa, - "OP ", - ) - elif physics_variables.ishape == 11: - po.ovarrf( - self.outfile, - "Elongation, X-point (calculated from aspect ratio via Menard 2016)", - "(kappa)", - physics_variables.kappa, - "OP ", - ) - else: - error_handling.idiags[0] = physics_variables.ishape - po.report_error(86) + # Issue #487. No idea what zeffai is. + # I haven't removed it as it is used in subroutine rether, + # (routine to find the equilibration power between the ions and electrons) + # po.ovarrf(self.outfile,'Mass weighted effective charge','(zeffai)',zeffai, 'OP ') - if physics_variables.ishape in [4, 5, 7]: - po.ovarrf( - self.outfile, - "Elongation, 95% surface (input value used)", - "(kappa95)", - physics_variables.kappa95, - "IP ", - ) - else: - po.ovarrf( - self.outfile, - "Elongation, 95% surface (calculated from kappa)", - "(kappa95)", - physics_variables.kappa95, - "OP ", - ) + po.ovarrf( + self.outfile, "Density profile factor", "(alphan)", physics_variables.alphan + ) + po.ovarin( + self.outfile, + "Plasma profile model", + "(ipedestal)", + physics_variables.ipedestal, + ) + if physics_variables.ipedestal >= 1: + if physics_variables.ne0 < physics_variables.neped: + error_handling.report_error(213) + + po.ocmmnt(self.outfile, "Pedestal profiles are used.") po.ovarrf( self.outfile, - "Elongation, area ratio calc.", - "(kappaa)", - physics_variables.kappaa, - "OP ", + "Density pedestal r/a location", + "(rhopedn)", + physics_variables.rhopedn, ) - - if physics_variables.ishape in [0, 2, 6, 8, 9, 10, 11]: - po.ovarrf( - self.outfile, - "Triangularity, X-point (input value used)", - "(triang)", - physics_variables.triang, - "IP ", - ) - elif physics_variables.ishape == 1: - po.ovarrf( + if physics_variables.fgwped >= 0e0: + po.ovarre( self.outfile, - "Triangularity, X-point (TART scaling)", - "(triang)", - physics_variables.triang, + "Electron density pedestal height (/m3)", + "(neped)", + physics_variables.neped, "OP ", ) else: - po.ovarrf( + po.ovarre( self.outfile, - "Triangularity, X-point (calculated from triang95)", - "(triang)", - physics_variables.triang, - "OP ", + "Electron density pedestal height (/m3)", + "(neped)", + physics_variables.neped, ) - if physics_variables.ishape in [3, 4, 5, 7]: - po.ovarrf( - self.outfile, - "Triangularity, 95% surface (input value used)", - "(triang95)", - physics_variables.triang95, - "IP ", + # This code is ODD# Don't change it# No explanation why fgwped and physics_variables.fgwsep + # must be assigned to their exisiting values# + fgwped_out = physics_variables.neped / physics_variables.dlimit[6] + fgwsep_out = physics_variables.nesep / physics_variables.dlimit[6] + if physics_variables.fgwped >= 0e0: + physics_variables.fgwped = ( + physics_variables.neped / physics_variables.dlimit[6] ) - else: - po.ovarrf( - self.outfile, - "Triangularity, 95% surface (calculated from triang)", - "(triang95)", - physics_variables.triang95, - "OP ", + if physics_variables.fgwsep >= 0e0: + physics_variables.fgwsep = ( + physics_variables.nesep / physics_variables.dlimit[6] ) - po.ovarrf( + po.ovarre( self.outfile, - "Plasma poloidal perimeter (m)", - "(pperim)", - physics_variables.pperim, - "OP ", + "Electron density at pedestal / nGW", + "(fgwped_out)", + fgwped_out, ) - po.ovarrf( self.outfile, - "Plasma cross-sectional area (m2)", - "(xarea)", - physics_variables.xarea, - "OP ", - ) - po.ovarre( - self.outfile, - "Plasma surface area (m2)", - "(sarea)", - physics_variables.sarea, - "OP ", + "Temperature pedestal r/a location", + "(rhopedt)", + physics_variables.rhopedt, ) - po.ovarre( + # Issue #413 Pedestal scaling + po.ovarin( self.outfile, - "Plasma volume (m3)", - "(vol)", - physics_variables.vol, - "OP ", + "Pedestal scaling switch", + "(ieped)", + physics_variables.ieped, ) + if physics_variables.ieped == 1: + po.ocmmnt( + self.outfile, + "Saarelma 6-parameter pedestal temperature scaling is ON", + ) - po.osubhd(self.outfile, "Current and Field :") - - if stellarator_variables.istell == 0: - if physics_variables.iprofile == 0: + if self.eped_warning() != "": po.ocmmnt( self.outfile, - "Consistency between q0,q,alphaj,rli,dnbeta is not enforced", + "WARNING: Pedestal parameters are outside the range of applicability of the scaling:", ) - else: po.ocmmnt( self.outfile, - "Consistency between q0,q,alphaj,rli,dnbeta is enforced", + "triang: 0.4 - 0.6; physics_variables.kappa: 1.5 - 2.0; plascur: 10 - 20 MA, physics_variables.rmajor: 7 - 11 m;", ) - - po.oblnkl(self.outfile) - po.ovarin( - self.outfile, - "Plasma current scaling law used", - "(icurr)", - physics_variables.icurr, - ) - - po.ovarrf( - self.outfile, - "Plasma current (MA)", - "(plascur/1D6)", - physics_variables.plascur / 1.0e6, - "OP ", - ) - - if physics_variables.iprofile == 1: - po.ovarrf( + po.ocmmnt( self.outfile, - "Current density profile factor", - "(alphaj)", - physics_variables.alphaj, - "OP ", + "rminor: 2 - 3.5 m; tesep: 0 - 0.5 keV; normalised_total_beta: 2 - 3; ", ) - else: - po.ovarrf( - self.outfile, - "Current density profile factor", - "(alphaj)", - physics_variables.alphaj, + print( + "WARNING: Pedestal parameters are outside the range of applicability of the scaling:" + ) + print( + "triang: 0.4 - 0.6; physics_variables.kappa: 1.5 - 2.0; plascur: 10 - 20 MA, physics_variables.rmajor: 7 - 11 m;" + ) + print( + "rminor: 2 - 3.5 m; tesep: 0 - 0.5 keV; normalised_total_beta: 2 - 3" ) + print(self.eped_warning()) + po.ovarrf( + self.outfile, + "Electron temp. pedestal height (keV)", + "(teped)", + physics_variables.teped, + ) + if any(numerics.icc == 78): po.ovarrf( self.outfile, - "Plasma internal inductance, li", - "(rli)", - physics_variables.rli, + "Electron temp. at separatrix (keV)", + "(tesep)", + physics_variables.tesep, "OP ", ) + else: po.ovarrf( self.outfile, - "Vertical field at plasma (T)", - "(bvert)", - physics_variables.bvert, - "OP ", + "Electron temp. at separatrix (keV)", + "(tesep)", + physics_variables.tesep, ) - po.ovarrf( + po.ovarre( self.outfile, - "Vacuum toroidal field at R (T)", - "(bt)", - physics_variables.bt, + "Electron density at separatrix (/m3)", + "(nesep)", + physics_variables.nesep, ) - po.ovarrf( + po.ovarre( self.outfile, - "Average poloidal field (T)", - "(bp)", - physics_variables.bp, - "OP ", + "Electron density at separatrix / nGW", + "(fgwsep_out)", + fgwsep_out, ) - po.ovarrf( + # Issue 558 - addition of constraint 76 to limit the value of nesep, in proportion with the ballooning parameter and Greenwald density + if any(numerics.icc == 76): + po.ovarre( self.outfile, - "Total field (sqrt(bp^2 + bt^2)) (T)", - "(btot)", - physics_variables.btot, - "OP ", + "Critical ballooning parameter value", + "(alpha_crit)", + physics_variables.alpha_crit, ) - - if stellarator_variables.istell == 0: - po.ovarrf( - self.outfile, "Safety factor on axis", "(q0)", physics_variables.q0 + po.ovarre( + self.outfile, + "Critical electron density at separatrix (/m3)", + "(nesep_crit)", + physics_variables.nesep_crit, ) - if physics_variables.icurr == 2: - po.ovarrf( - self.outfile, "Mean edge safety factor", "(q)", physics_variables.q - ) + po.ovarrf( + self.outfile, + "Temperature profile index", + "(alphat)", + physics_variables.alphat, + ) + po.ovarrf( + self.outfile, + "Temperature profile index beta", + "(tbeta)", + physics_variables.tbeta, + ) - po.ovarrf( + if stellarator_variables.istell == 0: + po.osubhd(self.outfile, "Density Limit using different models :") + po.ovarre( self.outfile, - "Safety factor at 95% flux surface", - "(q95)", - physics_variables.q95, + "Old ASDEX model", + "(dlimit(1))", + physics_variables.dlimit[0], + "OP ", ) - - po.ovarrf( + po.ovarre( self.outfile, - "Cylindrical safety factor (qcyl)", - "(qstar)", - physics_variables.qstar, + "Borrass ITER model I", + "(dlimit(2))", + physics_variables.dlimit[1], "OP ", ) - - if physics_variables.ishape == 1: - po.ovarrf( - self.outfile, - "Lower limit for edge safety factor q", - "(qlim)", - physics_variables.qlim, - "OP ", - ) - - else: - po.ovarrf( + po.ovarre( self.outfile, - "Rotational transform", - "(iotabar)", - stellarator_variables.iotabar, + "Borrass ITER model II", + "(dlimit(3))", + physics_variables.dlimit[2], + "OP ", + ) + po.ovarre( + self.outfile, + "JET edge radiation model", + "(dlimit(4))", + physics_variables.dlimit[3], + "OP ", + ) + po.ovarre( + self.outfile, + "JET simplified model", + "(dlimit(5))", + physics_variables.dlimit[4], + "OP ", + ) + po.ovarre( + self.outfile, + "Hugill-Murakami Mq model", + "(dlimit(6))", + physics_variables.dlimit[5], + "OP ", + ) + po.ovarre( + self.outfile, + "Greenwald model", + "(dlimit(7))", + physics_variables.dlimit[6], + "OP ", ) - po.osubhd(self.outfile, "Beta Information :") - - betath = ( - physics_variables.beta - physics_variables.betaft - physics_variables.betanb + po.osubhd(self.outfile, "Fuel Constituents :") + po.ovarrf( + self.outfile, "Deuterium fuel fraction", "(fdeut)", physics_variables.fdeut ) - gammaft = (physics_variables.betaft + physics_variables.betanb) / betath + po.ovarrf( + self.outfile, "Tritium fuel fraction", "(ftrit)", physics_variables.ftrit + ) + if physics_variables.fhe3 > 1.0e-3: + po.ovarrf( + self.outfile, "3-Helium fuel fraction", "(fhe3)", physics_variables.fhe3 + ) - po.ovarre(self.outfile, "Total plasma beta", "(beta)", physics_variables.beta) + po.osubhd(self.outfile, "Fusion Power :") po.ovarre( self.outfile, - "Total poloidal beta", - "(betap)", - physics_variables.betap, + "Total fusion power (MW)", + "(powfmw)", + physics_variables.powfmw, "OP ", ) po.ovarre( self.outfile, - "Total toroidal beta", - " ", - physics_variables.beta - * (physics_variables.btot / physics_variables.bt) ** 2, + " = D-T fusion power (MW)", + "(pdt)", + physics_variables.pdt, "OP ", ) po.ovarre( - self.outfile, "Fast alpha beta", "(betaft)", physics_variables.betaft, "OP " + self.outfile, + " + D-D fusion power (MW)", + "(pdd)", + physics_variables.pdd, + "OP ", ) po.ovarre( - self.outfile, "Beam ion beta", "(betanb)", physics_variables.betanb, "OP " + self.outfile, + " + D-He3 fusion power (MW)", + "(pdhe3)", + physics_variables.pdhe3, + "OP ", ) po.ovarre( self.outfile, - "(Fast alpha + beam physics_variables.beta)/(thermal physics_variables.beta)", - "(gammaft)", - gammaft, + "Alpha power: total (MW)", + "(palpmw)", + physics_variables.palpmw, "OP ", ) - - po.ovarre(self.outfile, "Thermal beta", " ", betath, "OP ") po.ovarre( self.outfile, - "Thermal poloidal beta", - " ", - betath * (physics_variables.btot / physics_variables.bp) ** 2, + "Alpha power: beam-plasma (MW)", + "(palpnb)", + physics_variables.palpnb, "OP ", ) po.ovarre( self.outfile, - "Thermal toroidal physics_variables.beta (= beta-exp)", - " ", - betath * (physics_variables.btot / physics_variables.bt) ** 2, + "Neutron power (MW)", + "(pneutmw)", + physics_variables.pneutmw, + "OP ", + ) + po.ovarre( + self.outfile, + "Charged particle power (excluding alphas) (MW)", + "(pchargemw)", + physics_variables.pchargemw, + "OP ", + ) + tot_power_plasma = ( + physics_variables.falpha * physics_variables.palpmw + + physics_variables.pchargemw + + physics_variables.pohmmw + + current_drive_variables.pinjmw + ) + po.ovarre( + self.outfile, + "Total power deposited in plasma (MW)", + "(tot_power_plasma)", + tot_power_plasma, "OP ", ) + # po.ovarre(self.outfile,'Total power deposited in plasma (MW)','()',falpha*palpmw+pchargemw+pohmmw+pinjmw, 'OP ') - po.ovarrf( + po.osubhd(self.outfile, "Radiation Power (excluding SOL):") + po.ovarre( self.outfile, - "2nd stability physics_variables.beta : beta_p / (R/a)", - "(eps*betap)", - physics_variables.eps * physics_variables.betap, + "Bremsstrahlung radiation power (MW)", + "(pbrempv*vol)", + physics_variables.pbrempv * physics_variables.vol, + "OP ", + ) + po.ovarre( + self.outfile, + "Line radiation power (MW)", + "(plinepv*vol)", + physics_variables.plinepv * physics_variables.vol, + "OP ", + ) + po.ovarre( + self.outfile, + "Synchrotron radiation power (MW)", + "(psyncpv*vol)", + physics_variables.psyncpv * physics_variables.vol, "OP ", ) po.ovarrf( self.outfile, - "2nd stability physics_variables.beta upper limit", - "(epbetmax)", - physics_variables.epbetmax, + "Synchrotron wall reflectivity factor", + "(ssync)", + physics_variables.ssync, ) - - if stellarator_variables.istell == 0: - if physics_variables.iprofile == 1: - po.ovarrf( - self.outfile, - "Beta g coefficient", - "(dnbeta)", - physics_variables.dnbeta, - "OP ", - ) - else: - po.ovarrf( - self.outfile, - "Beta g coefficient", - "(dnbeta)", - physics_variables.dnbeta, - ) - - po.ovarrf( - self.outfile, - "Normalised thermal beta", - " ", - 1.0e8 - * betath - * physics_variables.rminor - * physics_variables.bt - / physics_variables.plascur, - "OP ", - ) - - po.ovarrf( - self.outfile, - "Normalised total beta", - " ", - physics_variables.normalised_total_beta, - "OP ", - ) - - normalised_toroidal_beta = ( - physics_variables.normalised_total_beta - * (physics_variables.btot / physics_variables.bt) ** 2 - ) - po.ovarrf( - self.outfile, - "Normalised toroidal beta", - "(normalised_toroidal_beta)", - normalised_toroidal_beta, - "OP ", - ) - - if physics_variables.iculbl == 0: - po.ovarrf( - self.outfile, - "Limit on total beta", - "(betalim)", - physics_variables.betalim, - "OP ", - ) - elif physics_variables.iculbl == 1: - po.ovarrf( - self.outfile, - "Limit on thermal beta", - "(betalim)", - physics_variables.betalim, - "OP ", - ) - else: - po.ovarrf( - self.outfile, - "Limit on thermal + NB beta", - "(betalim)", - physics_variables.betalim, - "OP ", - ) - po.ovarre( self.outfile, - "Plasma thermal energy (J)", - " ", - 1.5e0 - * betath - * physics_variables.btot - * physics_variables.btot - / (2.0e0 * constants.rmu0) - * physics_variables.vol, - "OP ", + "Normalised minor radius defining 'core'", + "(coreradius)", + impurity_radiation_module.coreradius, ) - po.ovarre( self.outfile, - "Total plasma internal energy (J)", - "(total_plasma_internal_energy)", - physics_module.total_plasma_internal_energy, - "OP ", - ) - - po.osubhd(self.outfile, "Temperature and Density (volume averaged) :") - po.ovarrf( - self.outfile, "Electron temperature (keV)", "(te)", physics_variables.te - ) - po.ovarrf( - self.outfile, - "Electron temperature on axis (keV)", - "(te0)", - physics_variables.te0, - "OP ", - ) - po.ovarrf(self.outfile, "Ion temperature (keV)", "(ti)", physics_variables.ti) - po.ovarrf( - self.outfile, - "Ion temperature on axis (keV)", - "(ti0)", - physics_variables.ti0, - "OP ", - ) - po.ovarrf( - self.outfile, - "Electron temp., density weighted (keV)", - "(ten)", - physics_variables.ten, - "OP ", - ) - po.ovarre( - self.outfile, "Electron density (/m3)", "(dene)", physics_variables.dene + "Fraction of core radiation subtracted from P_L", + "(coreradiationfraction)", + impurity_radiation_module.coreradiationfraction, ) po.ovarre( self.outfile, - "Electron density on axis (/m3)", - "(ne0)", - physics_variables.ne0, + "Radiation power from inner zone (MW)", + "(pinnerzoneradmw)", + physics_variables.pinnerzoneradmw, "OP ", ) po.ovarre( self.outfile, - "Line-averaged electron density (/m3)", - "(dnla)", - physics_variables.dnla, + "Radiation power from outer zone (MW)", + "(pouterzoneradmw)", + physics_variables.pouterzoneradmw, "OP ", ) - if stellarator_variables.istell == 0: + if stellarator_variables.istell != 0: po.ovarre( self.outfile, - "Line-averaged electron density / Greenwald density", - "(dnla_gw)", - physics_variables.dnla / physics_variables.dlimit[6], + "SOL radiation power as imposed by f_rad (MW)", + "(psolradmw)", + physics_variables.psolradmw, "OP ", ) po.ovarre( self.outfile, - "Ion density (/m3)", - "(dnitot)", - physics_variables.dnitot, + "Total radiation power from inside LCFS (MW)", + "(pradmw)", + physics_variables.pradmw, "OP ", ) po.ovarre( - self.outfile, "Fuel density (/m3)", "(deni)", physics_variables.deni, "OP " + self.outfile, + "LCFS radiation fraction = total radiation in LCFS / total power deposited in plasma", + "(rad_fraction_LCFS)", + physics_module.rad_fraction_lcfs, + "OP ", ) po.ovarre( self.outfile, - "Total impurity density with Z > 2 (no He) (/m3)", - "(dnz)", - physics_variables.dnz, + "Nominal mean radiation load on inside surface of reactor (MW/m2)", + "(photon_wall)", + physics_variables.photon_wall, "OP ", ) po.ovarre( self.outfile, - "Helium ion density (thermalised ions only) (/m3)", - "(dnalp)", - physics_variables.dnalp, - "OP ", + "Peaking factor for radiation wall load", + "(peakfactrad)", + constraint_variables.peakfactrad, + "IP ", ) po.ovarre( self.outfile, - "Proton density (/m3)", - "(dnprot)", - physics_variables.dnprot, + "Maximum permitted radiation wall load (MW/m^2)", + "(maxradwallload)", + constraint_variables.maxradwallload, + "IP ", + ) + po.ovarre( + self.outfile, + "Peak radiation wall load (MW/m^2)", + "(peakradwallload)", + constraint_variables.peakradwallload, "OP ", ) - if physics_variables.protium > 1.0e-10: - po.ovarre( - self.outfile, - "Seeded protium density / electron density", - "(protium)", - physics_variables.protium, - ) - po.ovarre( self.outfile, - "Hot beam density (/m3)", - "(dnbeam)", - physics_variables.dnbeam, + "Fast alpha particle power incident on the first wall (MW)", + "(palpfwmw)", + physics_variables.palpfwmw, "OP ", ) po.ovarre( self.outfile, - "Density limit from scaling (/m3)", - "(dnelimt)", - physics_variables.dnelimt, + "Nominal mean neutron load on inside surface of reactor (MW/m2)", + "(wallmw)", + physics_variables.wallmw, "OP ", ) - if (numerics.ioptimz > 0) and (numerics.active_constraints[4]): + + if stellarator_variables.istell == 0: + po.oblnkl(self.outfile) po.ovarre( self.outfile, - "Density limit (enforced) (/m3)", - "(boundu(9)*dnelimt)", - numerics.boundu[8] * physics_variables.dnelimt, + "Power incident on the divertor targets (MW)", + "(ptarmw)", + physics_module.ptarmw, "OP ", ) - - po.ovarre( - self.outfile, - "Helium ion density (thermalised ions only) / electron density", - "(ralpne)", - physics_variables.ralpne, - ) - po.oblnkl(self.outfile) - - po.ocmmnt(self.outfile, "Impurities") - po.oblnkl(self.outfile) - po.ocmmnt(self.outfile, "Plasma ion densities / electron density:") - - for imp in range(impurity_radiation_module.nimp): - # MDK Update fimp, as this will make the ITV output work correctly. - impurity_radiation_module.fimp[ - imp - ] = impurity_radiation_module.impurity_arr_frac[imp] - str1 = ( - f2py_compatible_to_string( - impurity_radiation_module.impurity_arr_label[imp] - ) - + " concentration" + po.ovarre( + self.outfile, + "Fraction of power to the lower divertor", + "(ftar)", + physics_variables.ftar, + "IP ", ) - str2 = f"(fimp({imp+1:02}))" - # MDK Add output flag for H which is calculated. - if imp == 0: + po.ovarre( + self.outfile, + "Outboard side heat flux decay length (m)", + "(lambdaio)", + physics_module.lambdaio, + "OP ", + ) + if physics_variables.idivrt == 2: po.ovarre( - self.outfile, str1, str2, impurity_radiation_module.fimp[imp], "OP " + self.outfile, + "Midplane seperation of the two magnetic closed flux surfaces (m)", + "(drsep)", + physics_module.drsep, + "OP ", ) - else: - po.ovarre(self.outfile, str1, str2, impurity_radiation_module.fimp[imp]) - po.ovarre( - self.outfile, - "Average mass of all ions (amu)", - "(aion)", - physics_variables.aion, - "OP ", - ) - # MDK Say which impurity is varied, if iteration variable fimpvar (102) is turned on - # if (any(ixc == 102)) : - # call ovarst(self.outfile,'Impurity used as an iteration variable' , '', '"' // impurity_arr(impvar)%label // '"') - # po.ovarre(self.outfile,'Fractional density of variable impurity (ion / electron density)','(fimpvar)',fimpvar) - # - po.oblnkl(self.outfile) - po.ovarrf( - self.outfile, "Effective charge", "(zeff)", physics_variables.zeff, "OP " - ) - - # Issue #487. No idea what zeffai is. - # I haven't removed it as it is used in subroutine rether, - # (routine to find the equilibration power between the ions and electrons) - # po.ovarrf(self.outfile,'Mass weighted effective charge','(zeffai)',zeffai, 'OP ') - - po.ovarrf( - self.outfile, "Density profile factor", "(alphan)", physics_variables.alphan - ) - po.ovarin( - self.outfile, - "Plasma profile model", - "(ipedestal)", - physics_variables.ipedestal, - ) - - if physics_variables.ipedestal >= 1: - if physics_variables.ne0 < physics_variables.neped: - error_handling.report_error(213) - - po.ocmmnt(self.outfile, "Pedestal profiles are used.") - po.ovarrf( + po.ovarre( self.outfile, - "Density pedestal r/a location", - "(rhopedn)", - physics_variables.rhopedn, + "Fraction of power on the inner targets", + "(fio)", + physics_module.fio, + "OP ", ) - if physics_variables.fgwped >= 0e0: + po.ovarre( + self.outfile, + "Fraction of power incident on the lower inner target", + "(fLI)", + physics_module.fli, + "OP ", + ) + po.ovarre( + self.outfile, + "Fraction of power incident on the lower outer target", + "(fLO)", + physics_module.flo, + "OP ", + ) + if physics_variables.idivrt == 2: po.ovarre( self.outfile, - "Electron density pedestal height (/m3)", - "(neped)", - physics_variables.neped, + "Fraction of power incident on the upper inner target", + "(fUI)", + physics_module.fui, "OP ", ) - else: po.ovarre( self.outfile, - "Electron density pedestal height (/m3)", - "(neped)", - physics_variables.neped, - ) - - # This code is ODD# Don't change it# No explanation why fgwped and physics_variables.fgwsep - # must be assigned to their exisiting values# - fgwped_out = physics_variables.neped / physics_variables.dlimit[6] - fgwsep_out = physics_variables.nesep / physics_variables.dlimit[6] - if physics_variables.fgwped >= 0e0: - physics_variables.fgwped = ( - physics_variables.neped / physics_variables.dlimit[6] - ) - if physics_variables.fgwsep >= 0e0: - physics_variables.fgwsep = ( - physics_variables.nesep / physics_variables.dlimit[6] + "Fraction of power incident on the upper outer target", + "(fUO)", + physics_module.fuo, + "OP ", ) po.ovarre( self.outfile, - "Electron density at pedestal / nGW", - "(fgwped_out)", - fgwped_out, - ) - po.ovarrf( - self.outfile, - "Temperature pedestal r/a location", - "(rhopedt)", - physics_variables.rhopedt, - ) - # Issue #413 Pedestal scaling - po.ovarin( - self.outfile, - "Pedestal scaling switch", - "(ieped)", - physics_variables.ieped, + "Power incident on the lower inner target (MW)", + "(pLImw)", + physics_module.plimw, + "OP ", ) - if physics_variables.ieped == 1: - po.ocmmnt( - self.outfile, - "Saarelma 6-parameter pedestal temperature scaling is ON", - ) - - if self.eped_warning() != "": - po.ocmmnt( - self.outfile, - "WARNING: Pedestal parameters are outside the range of applicability of the scaling:", - ) - po.ocmmnt( - self.outfile, - "triang: 0.4 - 0.6; physics_variables.kappa: 1.5 - 2.0; plascur: 10 - 20 MA, physics_variables.rmajor: 7 - 11 m;", - ) - po.ocmmnt( - self.outfile, - "rminor: 2 - 3.5 m; tesep: 0 - 0.5 keV; normalised_total_beta: 2 - 3; ", - ) - print( - "WARNING: Pedestal parameters are outside the range of applicability of the scaling:" - ) - print( - "triang: 0.4 - 0.6; physics_variables.kappa: 1.5 - 2.0; plascur: 10 - 20 MA, physics_variables.rmajor: 7 - 11 m;" - ) - print( - "rminor: 2 - 3.5 m; tesep: 0 - 0.5 keV; normalised_total_beta: 2 - 3" - ) - print(self.eped_warning()) - - po.ovarrf( + po.ovarre( self.outfile, - "Electron temp. pedestal height (keV)", - "(teped)", - physics_variables.teped, + "Power incident on the lower outer target (MW)", + "(pLOmw)", + physics_module.plomw, + "OP ", ) - if any(numerics.icc == 78): - po.ovarrf( + if physics_variables.idivrt == 2: + po.ovarre( self.outfile, - "Electron temp. at separatrix (keV)", - "(tesep)", - physics_variables.tesep, + "Power incident on the upper innner target (MW)", + "(pUImw)", + physics_module.puimw, "OP ", ) - else: - po.ovarrf( + po.ovarre( self.outfile, - "Electron temp. at separatrix (keV)", - "(tesep)", - physics_variables.tesep, + "Power incident on the upper outer target (MW)", + "(pUOmw)", + physics_module.puomw, + "OP ", ) - po.ovarre( - self.outfile, - "Electron density at separatrix (/m3)", - "(nesep)", - physics_variables.nesep, - ) - po.ovarre( - self.outfile, - "Electron density at separatrix / nGW", - "(fgwsep_out)", - fgwsep_out, - ) - - # Issue 558 - addition of constraint 76 to limit the value of nesep, in proportion with the ballooning parameter and Greenwald density - if any(numerics.icc == 76): - po.ovarre( - self.outfile, - "Critical ballooning parameter value", - "(alpha_crit)", - physics_variables.alpha_crit, - ) - po.ovarre( - self.outfile, - "Critical electron density at separatrix (/m3)", - "(nesep_crit)", - physics_variables.nesep_crit, - ) - + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Ohmic heating power (MW)", + "(pohmmw)", + physics_variables.pohmmw, + "OP ", + ) po.ovarrf( self.outfile, - "Temperature profile index", - "(alphat)", - physics_variables.alphat, + "Fraction of alpha power deposited in plasma", + "(falpha)", + physics_variables.falpha, + "OP ", ) po.ovarrf( self.outfile, - "Temperature profile index beta", - "(tbeta)", - physics_variables.tbeta, + "Fraction of alpha power to electrons", + "(falpe)", + physics_variables.falpe, + "OP ", + ) + po.ovarrf( + self.outfile, + "Fraction of alpha power to ions", + "(falpi)", + physics_variables.falpi, + "OP ", + ) + po.ovarre( + self.outfile, + "Ion transport (MW)", + "(ptrimw)", + physics_variables.ptrimw, + "OP ", + ) + po.ovarre( + self.outfile, + "Electron transport (MW)", + "(ptremw)", + physics_variables.ptremw, + "OP ", + ) + po.ovarre( + self.outfile, + "Injection power to ions (MW)", + "(pinjimw)", + current_drive_variables.pinjimw, + "OP ", + ) + po.ovarre( + self.outfile, + "Injection power to electrons (MW)", + "(pinjemw)", + current_drive_variables.pinjemw, + "OP ", ) + if physics_variables.ignite == 1: + po.ocmmnt(self.outfile, " (Injected power only used for start-up phase)") - if stellarator_variables.istell == 0: - po.osubhd(self.outfile, "Density Limit using different models :") - po.ovarre( - self.outfile, - "Old ASDEX model", - "(dlimit(1))", - physics_variables.dlimit[0], - "OP ", - ) - po.ovarre( - self.outfile, - "Borrass ITER model I", - "(dlimit(2))", - physics_variables.dlimit[1], + po.ovarin( + self.outfile, + "Ignited plasma switch (0=not ignited, 1=ignited)", + "(ignite)", + physics_variables.ignite, + ) + + po.oblnkl(self.outfile) + po.ovarre( + self.outfile, + "Power into divertor zone via charged particles (MW)", + "(pdivt)", + physics_variables.pdivt, + "OP ", + ) + + if physics_variables.pdivt <= 0.001e0: + error_handling.fdiags[0] = physics_variables.pdivt + error_handling.report_error(87) + po.oblnkl(self.outfile) + po.ocmmnt( + self.outfile, " BEWARE: possible problem with high radiation power" + ) + po.ocmmnt( + self.outfile, " Power into divertor zone is unrealistic;" + ) + po.ocmmnt(self.outfile, " divertor calculations will be nonsense#") + po.ocmmnt( + self.outfile, " Set constraint 17 (Radiation fraction upper limit)." + ) + po.oblnkl(self.outfile) + + if physics_variables.idivrt == 2: + # Double null divertor configuration + po.ovarre( + self.outfile, + "Pdivt / R ratio (MW/m) (On peak divertor)", + "(pdivmax/physics_variables.rmajor)", + physics_variables.pdivmax / physics_variables.rmajor, "OP ", ) po.ovarre( self.outfile, - "Borrass ITER model II", - "(dlimit(3))", - physics_variables.dlimit[2], + "Pdivt Bt / qAR ratio (MWT/m) (On peak divertor)", + "(pdivmaxbt/qar)", + ( + (physics_variables.pdivmax * physics_variables.bt) + / ( + physics_variables.q95 + * physics_variables.aspect + * physics_variables.rmajor + ) + ), "OP ", ) + else: + # Single null divertor configuration po.ovarre( self.outfile, - "JET edge radiation model", - "(dlimit(4))", - physics_variables.dlimit[3], + "Psep / R ratio (MW/m)", + "(pdivt/rmajor)", + physics_variables.pdivt / physics_variables.rmajor, "OP ", ) po.ovarre( self.outfile, - "JET simplified model", - "(dlimit(5))", - physics_variables.dlimit[4], + "Psep Bt / qAR ratio (MWT/m)", + "(pdivtbt/qar)", + ( + (physics_variables.pdivt * physics_variables.bt) + / ( + physics_variables.q95 + * physics_variables.aspect + * physics_variables.rmajor + ) + ), "OP ", ) + + if stellarator_variables.istell == 0: + po.osubhd(self.outfile, "H-mode Power Threshold Scalings :") + po.ovarre( self.outfile, - "Hugill-Murakami Mq model", - "(dlimit(6))", - physics_variables.dlimit[5], + "ITER 1996 scaling: nominal (MW)", + "(pthrmw(1))", + physics_variables.pthrmw[0], "OP ", ) po.ovarre( self.outfile, - "Greenwald model", - "(dlimit(7))", - physics_variables.dlimit[6], + "ITER 1996 scaling: upper bound (MW)", + "(pthrmw(2))", + physics_variables.pthrmw[1], "OP ", ) - - po.osubhd(self.outfile, "Fuel Constituents :") - po.ovarrf( - self.outfile, "Deuterium fuel fraction", "(fdeut)", physics_variables.fdeut - ) - po.ovarrf( - self.outfile, "Tritium fuel fraction", "(ftrit)", physics_variables.ftrit - ) - if physics_variables.fhe3 > 1.0e-3: - po.ovarrf( - self.outfile, "3-Helium fuel fraction", "(fhe3)", physics_variables.fhe3 + po.ovarre( + self.outfile, + "ITER 1996 scaling: lower bound (MW)", + "(pthrmw(3))", + physics_variables.pthrmw[2], + "OP ", ) - - po.osubhd(self.outfile, "Fusion Power :") - po.ovarre( - self.outfile, - "Total fusion power (MW)", - "(powfmw)", - physics_variables.powfmw, - "OP ", - ) - po.ovarre( - self.outfile, - " = D-T fusion power (MW)", - "(pdt)", - physics_variables.pdt, - "OP ", - ) - po.ovarre( - self.outfile, - " + D-D fusion power (MW)", - "(pdd)", - physics_variables.pdd, - "OP ", - ) - po.ovarre( - self.outfile, - " + D-He3 fusion power (MW)", - "(pdhe3)", - physics_variables.pdhe3, - "OP ", - ) - po.ovarre( - self.outfile, - "Alpha power: total (MW)", - "(palpmw)", - physics_variables.palpmw, - "OP ", - ) - po.ovarre( - self.outfile, - "Alpha power: beam-plasma (MW)", - "(palpnb)", - physics_variables.palpnb, - "OP ", - ) - po.ovarre( - self.outfile, - "Neutron power (MW)", - "(pneutmw)", - physics_variables.pneutmw, - "OP ", - ) - po.ovarre( - self.outfile, - "Charged particle power (excluding alphas) (MW)", - "(pchargemw)", - physics_variables.pchargemw, - "OP ", - ) - tot_power_plasma = ( - physics_variables.falpha * physics_variables.palpmw - + physics_variables.pchargemw - + physics_variables.pohmmw - + current_drive_variables.pinjmw - ) - po.ovarre( - self.outfile, - "Total power deposited in plasma (MW)", - "(tot_power_plasma)", - tot_power_plasma, - "OP ", - ) - # po.ovarre(self.outfile,'Total power deposited in plasma (MW)','()',falpha*palpmw+pchargemw+pohmmw+pinjmw, 'OP ') - - po.osubhd(self.outfile, "Radiation Power (excluding SOL):") - po.ovarre( - self.outfile, - "Bremsstrahlung radiation power (MW)", - "(pbrempv*vol)", - physics_variables.pbrempv * physics_variables.vol, - "OP ", - ) - po.ovarre( - self.outfile, - "Line radiation power (MW)", - "(plinepv*vol)", - physics_variables.plinepv * physics_variables.vol, - "OP ", - ) - po.ovarre( - self.outfile, - "Synchrotron radiation power (MW)", - "(psyncpv*vol)", - physics_variables.psyncpv * physics_variables.vol, - "OP ", - ) - po.ovarrf( - self.outfile, - "Synchrotron wall reflectivity factor", - "(ssync)", - physics_variables.ssync, - ) - po.ovarre( - self.outfile, - "Normalised minor radius defining 'core'", - "(coreradius)", - impurity_radiation_module.coreradius, - ) - po.ovarre( - self.outfile, - "Fraction of core radiation subtracted from P_L", - "(coreradiationfraction)", - impurity_radiation_module.coreradiationfraction, - ) - po.ovarre( - self.outfile, - "Radiation power from inner zone (MW)", - "(pinnerzoneradmw)", - physics_variables.pinnerzoneradmw, - "OP ", - ) - po.ovarre( - self.outfile, - "Radiation power from outer zone (MW)", - "(pouterzoneradmw)", - physics_variables.pouterzoneradmw, - "OP ", - ) - - if stellarator_variables.istell != 0: po.ovarre( self.outfile, - "SOL radiation power as imposed by f_rad (MW)", - "(psolradmw)", - physics_variables.psolradmw, + "ITER 1997 scaling (1) (MW)", + "(pthrmw(4))", + physics_variables.pthrmw[3], "OP ", ) - - po.ovarre( - self.outfile, - "Total radiation power from inside LCFS (MW)", - "(pradmw)", - physics_variables.pradmw, - "OP ", - ) - po.ovarre( - self.outfile, - "LCFS radiation fraction = total radiation in LCFS / total power deposited in plasma", - "(rad_fraction_LCFS)", - physics_module.rad_fraction_lcfs, - "OP ", - ) - po.ovarre( - self.outfile, - "Nominal mean radiation load on inside surface of reactor (MW/m2)", - "(photon_wall)", - physics_variables.photon_wall, - "OP ", - ) - po.ovarre( - self.outfile, - "Peaking factor for radiation wall load", - "(peakfactrad)", - constraint_variables.peakfactrad, - "IP ", - ) - po.ovarre( - self.outfile, - "Maximum permitted radiation wall load (MW/m^2)", - "(maxradwallload)", - constraint_variables.maxradwallload, - "IP ", - ) - po.ovarre( - self.outfile, - "Peak radiation wall load (MW/m^2)", - "(peakradwallload)", - constraint_variables.peakradwallload, - "OP ", - ) - po.ovarre( - self.outfile, - "Fast alpha particle power incident on the first wall (MW)", - "(palpfwmw)", - physics_variables.palpfwmw, - "OP ", - ) - po.ovarre( - self.outfile, - "Nominal mean neutron load on inside surface of reactor (MW/m2)", - "(wallmw)", - physics_variables.wallmw, - "OP ", - ) - - if stellarator_variables.istell == 0: - po.oblnkl(self.outfile) po.ovarre( self.outfile, - "Power incident on the divertor targets (MW)", - "(ptarmw)", - physics_module.ptarmw, + "ITER 1997 scaling (2) (MW)", + "(pthrmw(5))", + physics_variables.pthrmw[4], "OP ", ) po.ovarre( self.outfile, - "Fraction of power to the lower divertor", - "(ftar)", - physics_variables.ftar, - "IP ", + "Martin 2008 scaling: nominal (MW)", + "(pthrmw(6))", + physics_variables.pthrmw[5], + "OP ", ) po.ovarre( self.outfile, - "Outboard side heat flux decay length (m)", - "(lambdaio)", - physics_module.lambdaio, + "Martin 2008 scaling: 95% upper bound (MW)", + "(pthrmw(7))", + physics_variables.pthrmw[6], "OP ", ) - if physics_variables.idivrt == 2: - po.ovarre( - self.outfile, - "Midplane seperation of the two magnetic closed flux surfaces (m)", - "(drsep)", - physics_module.drsep, - "OP ", - ) - po.ovarre( self.outfile, - "Fraction of power on the inner targets", - "(fio)", - physics_module.fio, + "Martin 2008 scaling: 95% lower bound (MW)", + "(pthrmw(8))", + physics_variables.pthrmw[7], "OP ", ) po.ovarre( self.outfile, - "Fraction of power incident on the lower inner target", - "(fLI)", - physics_module.fli, + "Snipes 2000 scaling: nominal (MW)", + "(pthrmw(9))", + physics_variables.pthrmw[8], "OP ", ) po.ovarre( self.outfile, - "Fraction of power incident on the lower outer target", - "(fLO)", - physics_module.flo, + "Snipes 2000 scaling: upper bound (MW)", + "(pthrmw(10))", + physics_variables.pthrmw[9], "OP ", ) - if physics_variables.idivrt == 2: - po.ovarre( - self.outfile, - "Fraction of power incident on the upper inner target", - "(fUI)", - physics_module.fui, - "OP ", - ) - po.ovarre( - self.outfile, - "Fraction of power incident on the upper outer target", - "(fUO)", - physics_module.fuo, - "OP ", - ) - po.ovarre( self.outfile, - "Power incident on the lower inner target (MW)", - "(pLImw)", - physics_module.plimw, + "Snipes 2000 scaling: lower bound (MW)", + "(pthrmw(11))", + physics_variables.pthrmw[10], "OP ", ) po.ovarre( self.outfile, - "Power incident on the lower outer target (MW)", - "(pLOmw)", - physics_module.plomw, + "Snipes 2000 scaling (closed divertor): nominal (MW)", + "(pthrmw(12))", + physics_variables.pthrmw[11], "OP ", ) - if physics_variables.idivrt == 2: - po.ovarre( - self.outfile, - "Power incident on the upper innner target (MW)", - "(pUImw)", - physics_module.puimw, - "OP ", - ) - po.ovarre( - self.outfile, - "Power incident on the upper outer target (MW)", - "(pUOmw)", - physics_module.puomw, - "OP ", - ) - - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Ohmic heating power (MW)", - "(pohmmw)", - physics_variables.pohmmw, - "OP ", - ) - po.ovarrf( - self.outfile, - "Fraction of alpha power deposited in plasma", - "(falpha)", - physics_variables.falpha, - "OP ", - ) - po.ovarrf( - self.outfile, - "Fraction of alpha power to electrons", - "(falpe)", - physics_variables.falpe, - "OP ", - ) - po.ovarrf( - self.outfile, - "Fraction of alpha power to ions", - "(falpi)", - physics_variables.falpi, - "OP ", - ) - po.ovarre( - self.outfile, - "Ion transport (MW)", - "(ptrimw)", - physics_variables.ptrimw, - "OP ", - ) - po.ovarre( - self.outfile, - "Electron transport (MW)", - "(ptremw)", - physics_variables.ptremw, - "OP ", - ) - po.ovarre( - self.outfile, - "Injection power to ions (MW)", - "(pinjimw)", - current_drive_variables.pinjimw, - "OP ", - ) - po.ovarre( - self.outfile, - "Injection power to electrons (MW)", - "(pinjemw)", - current_drive_variables.pinjemw, - "OP ", - ) - if physics_variables.ignite == 1: - po.ocmmnt(self.outfile, " (Injected power only used for start-up phase)") - - po.ovarin( - self.outfile, - "Ignited plasma switch (0=not ignited, 1=ignited)", - "(ignite)", - physics_variables.ignite, - ) - - po.oblnkl(self.outfile) - po.ovarre( - self.outfile, - "Power into divertor zone via charged particles (MW)", - "(pdivt)", - physics_variables.pdivt, - "OP ", - ) - - if physics_variables.pdivt <= 0.001e0: - error_handling.fdiags[0] = physics_variables.pdivt - error_handling.report_error(87) - po.oblnkl(self.outfile) - po.ocmmnt( - self.outfile, " BEWARE: possible problem with high radiation power" - ) - po.ocmmnt( - self.outfile, " Power into divertor zone is unrealistic;" - ) - po.ocmmnt(self.outfile, " divertor calculations will be nonsense#") - po.ocmmnt( - self.outfile, " Set constraint 17 (Radiation fraction upper limit)." - ) - po.oblnkl(self.outfile) - - if physics_variables.idivrt == 2: - # Double null divertor configuration po.ovarre( self.outfile, - "Pdivt / R ratio (MW/m) (On peak divertor)", - "(pdivmax/physics_variables.rmajor)", - physics_variables.pdivmax / physics_variables.rmajor, + "Snipes 2000 scaling (closed divertor): upper bound (MW)", + "(pthrmw(13))", + physics_variables.pthrmw[12], "OP ", ) po.ovarre( self.outfile, - "Pdivt Bt / qAR ratio (MWT/m) (On peak divertor)", - "(pdivmaxbt/qar)", - ( - (physics_variables.pdivmax * physics_variables.bt) - / ( - physics_variables.q95 - * physics_variables.aspect - * physics_variables.rmajor - ) - ), - "OP ", - ) - else: - # Single null divertor configuration - po.ovarre( - self.outfile, - "Psep / R ratio (MW/m)", - "(pdivt/rmajor)", - physics_variables.pdivt / physics_variables.rmajor, - "OP ", - ) - po.ovarre( - self.outfile, - "Psep Bt / qAR ratio (MWT/m)", - "(pdivtbt/qar)", - ( - (physics_variables.pdivt * physics_variables.bt) - / ( - physics_variables.q95 - * physics_variables.aspect - * physics_variables.rmajor - ) - ), - "OP ", - ) - - if stellarator_variables.istell == 0: - po.osubhd(self.outfile, "H-mode Power Threshold Scalings :") - - po.ovarre( - self.outfile, - "ITER 1996 scaling: nominal (MW)", - "(pthrmw(1))", - physics_variables.pthrmw[0], - "OP ", - ) - po.ovarre( - self.outfile, - "ITER 1996 scaling: upper bound (MW)", - "(pthrmw(2))", - physics_variables.pthrmw[1], - "OP ", - ) - po.ovarre( - self.outfile, - "ITER 1996 scaling: lower bound (MW)", - "(pthrmw(3))", - physics_variables.pthrmw[2], - "OP ", - ) - po.ovarre( - self.outfile, - "ITER 1997 scaling (1) (MW)", - "(pthrmw(4))", - physics_variables.pthrmw[3], - "OP ", - ) - po.ovarre( - self.outfile, - "ITER 1997 scaling (2) (MW)", - "(pthrmw(5))", - physics_variables.pthrmw[4], - "OP ", - ) - po.ovarre( - self.outfile, - "Martin 2008 scaling: nominal (MW)", - "(pthrmw(6))", - physics_variables.pthrmw[5], - "OP ", - ) - po.ovarre( - self.outfile, - "Martin 2008 scaling: 95% upper bound (MW)", - "(pthrmw(7))", - physics_variables.pthrmw[6], - "OP ", - ) - po.ovarre( - self.outfile, - "Martin 2008 scaling: 95% lower bound (MW)", - "(pthrmw(8))", - physics_variables.pthrmw[7], - "OP ", - ) - po.ovarre( - self.outfile, - "Snipes 2000 scaling: nominal (MW)", - "(pthrmw(9))", - physics_variables.pthrmw[8], - "OP ", - ) - po.ovarre( - self.outfile, - "Snipes 2000 scaling: upper bound (MW)", - "(pthrmw(10))", - physics_variables.pthrmw[9], - "OP ", - ) - po.ovarre( - self.outfile, - "Snipes 2000 scaling: lower bound (MW)", - "(pthrmw(11))", - physics_variables.pthrmw[10], - "OP ", - ) - po.ovarre( - self.outfile, - "Snipes 2000 scaling (closed divertor): nominal (MW)", - "(pthrmw(12))", - physics_variables.pthrmw[11], - "OP ", - ) - po.ovarre( - self.outfile, - "Snipes 2000 scaling (closed divertor): upper bound (MW)", - "(pthrmw(13))", - physics_variables.pthrmw[12], - "OP ", - ) - po.ovarre( - self.outfile, - "Snipes 2000 scaling (closed divertor): lower bound (MW)", - "(pthrmw(14))", - physics_variables.pthrmw[13], + "Snipes 2000 scaling (closed divertor): lower bound (MW)", + "(pthrmw(14))", + physics_variables.pthrmw[13], "OP ", ) po.ovarre( @@ -4236,7 +3860,7 @@ def outplas(self): po.ocmmnt( self.outfile, - f"Confinement scaling law: {f2py_compatible_to_string(physics_variables.tauscl[physics_variables.isc-1])}", + f"Confinement scaling law: {f2py_compatible_to_string(physics_variables.tauscl[physics_variables.isc - 1])}", ) po.ovarrf( @@ -4524,202 +4148,543 @@ def outplas(self): if current_drive_variables.diacf_scene > 0.01e0: error_handling.report_error(244) - elif physics_variables.idia == 1: - po.ocmmnt( - self.outfile, " (Hender diamagnetic current fraction scaling used)" - ) - elif physics_variables.idia == 2: - po.ocmmnt( - self.outfile, " (SCENE diamagnetic current fraction scaling used)" - ) + elif physics_variables.idia == 1: + po.ocmmnt( + self.outfile, " (Hender diamagnetic current fraction scaling used)" + ) + elif physics_variables.idia == 2: + po.ocmmnt( + self.outfile, " (SCENE diamagnetic current fraction scaling used)" + ) + + if physics_variables.ips == 0: + po.ocmmnt( + self.outfile, " Pfirsch-Schluter current fraction not calculated" + ) + elif physics_variables.ips == 1: + po.ocmmnt( + self.outfile, + " (SCENE Pfirsch-Schluter current fraction scaling used)", + ) + + po.ovarrf( + self.outfile, + "Bootstrap fraction (enforced)", + "(bootipf.)", + current_drive_variables.bootipf, + "OP ", + ) + po.ovarrf( + self.outfile, + "Diamagnetic fraction (enforced)", + "(diaipf.)", + current_drive_variables.diaipf, + "OP ", + ) + po.ovarrf( + self.outfile, + "Pfirsch-Schlueter fraction (enforced)", + "(psipf.)", + current_drive_variables.psipf, + "OP ", + ) + + po.ovarre( + self.outfile, + "Loop voltage during burn (V)", + "(vburn)", + physics_variables.plascur + * physics_variables.rplas + * physics_variables.facoh, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma resistance (ohm)", + "(rplas)", + physics_variables.rplas, + "OP ", + ) + + po.ovarre( + self.outfile, + "Resistive diffusion time (s)", + "(res_time)", + physics_variables.res_time, + "OP ", + ) + po.ovarre( + self.outfile, + "Plasma inductance (H)", + "(rlp)", + physics_variables.rlp, + "OP ", + ) + po.ovarrf( + self.outfile, + "Coefficient for sawtooth effects on burn V-s requirement", + "(csawth)", + physics_variables.csawth, + ) + + po.osubhd(self.outfile, "Fuelling :") + po.ovarre( + self.outfile, + "Ratio of He and pellet particle confinement times", + "(tauratio)", + physics_variables.tauratio, + ) + po.ovarre( + self.outfile, + "Fuelling rate (nucleus-pairs/s)", + "(qfuel)", + physics_variables.qfuel, + "OP ", + ) + po.ovarre( + self.outfile, + "Fuel burn-up rate (reactions/s)", + "(rndfuel)", + physics_variables.rndfuel, + "OP ", + ) + po.ovarrf( + self.outfile, + "Burn-up fraction", + "(burnup)", + physics_variables.burnup, + "OP ", + ) + + if any(numerics.icc == 78): + po.osubhd(self.outfile, "Reinke Criterion :") + po.ovarin( + self.outfile, + "index of impurity to be iterated for divertor detachment", + "(impvardiv)", + reinke_variables.impvardiv, + ) + po.ovarre( + self.outfile, + "Minimum Impurity fraction from Reinke", + "(fzmin)", + reinke_variables.fzmin, + "OP ", + ) + po.ovarre( + self.outfile, + "Actual Impurity fraction", + "(fzactual)", + reinke_variables.fzactual, + ) + + def igmarcal(self): + """Routine to calculate ignition margin + author: P J Knight, CCFE, Culham Science Centre + outfile : input integer : Fortran output unit identifier + This routine calculates the ignition margin at the final point + with different scalings. + AEA FUS 251: A User's Guide to the PROCESS Systems Code + """ + + po.oheadr(self.outfile, "Energy confinement times, and required H-factors :") + po.ocmmnt( + self.outfile, + f"{'':>5}{'scaling law':<25}{'confinement time (s)':<25}H-factor for", + ) + po.ocmmnt( + self.outfile, + f"{'':>34}{'for H = 1':<20}power balance", + ) + + for iisc in range(32, 48): + ( + physics_variables.kappaa, + ptrez, + ptriz, + taueez, + taueiz, + taueffz, + powerhtz, + ) = self.pcond( + physics_variables.afuel, + physics_variables.palpmw, + physics_variables.aspect, + physics_variables.bt, + physics_variables.dnitot, + physics_variables.dene, + physics_variables.dnla, + physics_variables.eps, + 1.0, + physics_variables.iinvqd, + iisc, + physics_variables.ignite, + physics_variables.kappa, + physics_variables.kappa95, + physics_variables.pchargemw, + current_drive_variables.pinjmw, + physics_variables.plascur, + physics_variables.pcoreradpv, + physics_variables.rmajor, + physics_variables.rminor, + physics_variables.te, + physics_variables.ten, + physics_variables.tin, + physics_variables.q, + physics_variables.qstar, + physics_variables.vol, + physics_variables.xarea, + physics_variables.zeff, + ) + + physics_variables.hfac[iisc - 1] = self.fhfac(iisc) + + po.ocmmnt( + self.outfile, + f"{'':>2}{f2py_compatible_to_string(physics_variables.tauscl[iisc - 1]):<32}" + f"{taueez:<26.3f}{physics_variables.hfac[iisc - 1]:.3f}", + ) + + @staticmethod + def bootstrap_fraction_iter89( + aspect, beta, bt, cboot, plascur, q95, q0, rmajor, vol + ): + """Original ITER calculation of bootstrap-driven fraction + of the plasma current. + author: P J Knight, CCFE, Culham Science Centre + aspect : input real : plasma aspect ratio + beta : input real : plasma total beta + bt : input real : toroidal field on axis (T) + cboot : input real : bootstrap current fraction multiplier + plascur : input real : plasma current (A) + q95 : input real : safety factor at 95% surface + q0 : input real : central safety factor + rmajor : input real : plasma major radius (m) + vol : input real : plasma volume (m3) + This routine performs the original ITER calculation of the + plasma current bootstrap fraction. + ITER Physics Design Guidelines: 1989 [IPDG89], N. A. Uckan et al, + ITER Documentation Series No.10, IAEA/ITER/DS/10, IAEA, Vienna, 1990 + """ + xbs = min(10, q95 / q0) + cbs = cboot * (1.32 - 0.235 * xbs + 0.0185 * xbs**2) + bpbs = ( + constants.rmu0 + * plascur + / (2 * np.pi * np.sqrt(vol / (2 * np.pi**2 * rmajor))) + ) + betapbs = beta * bt**2 / bpbs**2 + + if betapbs <= 0.0: # only possible if beta <= 0.0 + return 0.0 + return cbs * (betapbs / np.sqrt(aspect)) ** 1.3 + + @staticmethod + def bootstrap_fraction_wilson( + alphaj, alphap, alphat, betpth, q0, qpsi, rmajor, rminor + ): + """Bootstrap current fraction from Wilson et al scaling + author: P J Knight, CCFE, Culham Science Centre + alphaj : input real : current profile index + alphap : input real : pressure profile index + alphat : input real : temperature profile index + beta : input real : total beta + betpth : input real : thermal component of poloidal beta + q0 : input real : safety factor on axis + qpsi : input real : edge safety factor + rmajor : input real : major radius (m) + rminor : input real : minor radius (m) + This function calculates the bootstrap current fraction + using the numerically fitted algorithm written by Howard Wilson. + AEA FUS 172: Physics Assessment for the European Reactor Study + H. R. Wilson, Nuclear Fusion 32 (1992) 257 + """ + term1 = np.log(0.5) + term2 = np.log(q0 / qpsi) + + termp = 1.0 - 0.5 ** (1.0 / alphap) + termt = 1.0 - 0.5 ** (1.0 / alphat) + termj = 1.0 - 0.5 ** (1.0 / alphaj) + + alfpnw = term1 / np.log(np.log((q0 + (qpsi - q0) * termp) / qpsi) / term2) + alftnw = term1 / np.log(np.log((q0 + (qpsi - q0) * termt) / qpsi) / term2) + aj = term1 / np.log(np.log((q0 + (qpsi - q0) * termj) / qpsi) / term2) + + # Crude check for NaN errors or other illegal values... + + if np.isnan(aj) or np.isnan(alfpnw) or np.isnan(alftnw) or aj < 0: + error_handling.fdiags[0] = aj + error_handling.fdiags[1] = alfpnw + error_handling.fdiags[2] = alftnw + error_handling.fdiags[3] = aj + + error_handling.report_error(76) + + # Ratio of ionic charge to electron charge + + z = 1.0 + + # Inverse aspect ratio: r2 = maximum plasma radius, r1 = minimum + + r2 = rmajor + rminor + r1 = rmajor - rminor + eps1 = (r2 - r1) / (r2 + r1) + + # Coefficients fitted using least squares techniques + saj = np.sqrt(aj) + + a = np.array( + [ + 1.41 * (1.0 - 0.28 * saj) * (1.0 + 0.12 / z), + 0.36 * (1.0 - 0.59 * saj) * (1.0 + 0.8 / z), + -0.27 * (1.0 - 0.47 * saj) * (1.0 + 3.0 / z), + 0.0053 * (1.0 + 5.0 / z), + -0.93 * (1.0 - 0.34 * saj) * (1.0 + 0.15 / z), + -0.26 * (1.0 - 0.57 * saj) * (1.0 - 0.27 * z), + 0.064 * (1.0 - 0.6 * aj + 0.15 * aj * aj) * (1.0 + 7.6 / z), + -0.0011 * (1.0 + 9.0 / z), + -0.33 * (1.0 - aj + 0.33 * aj * aj), + -0.26 * (1.0 - 0.87 / saj - 0.16 * aj), + -0.14 * (1.0 - 1.14 / saj - 0.45 * saj), + -0.0069, + ] + ) + + seps1 = np.sqrt(eps1) + + b = np.array( + [ + 1.0, + alfpnw, + alftnw, + alfpnw * alftnw, + seps1, + alfpnw * seps1, + alftnw * seps1, + alfpnw * alftnw * seps1, + eps1, + alfpnw * eps1, + alftnw * eps1, + alfpnw * alftnw * eps1, + ] + ) + + # Empirical bootstrap current fraction + return seps1 * betpth * (a * b).sum() + + @staticmethod + def bootstrap_fraction_nevins( + alphan, + alphat, + betat, + bt, + dene, + plascur, + q95, + q0, + rmajor, + rminor, + ten, + zeff, + ): + """Bootstrap current fraction from Nevins et al scaling + author: P J Knight, CCFE, Culham Science Centre + alphan : input real : density profile index + alphat : input real : temperature profile index + betat : input real : total plasma beta (with respect to the toroidal + field) + bt : input real : toroidal field on axis (T) + dene : input real : electron density (/m3) + plascur: input real : plasma current (A) + q0 : input real : central safety factor + q95 : input real : safety factor at 95% surface + rmajor : input real : plasma major radius (m) + rminor : input real : plasma minor radius (m) + ten : input real : density weighted average plasma temperature (keV) + zeff : input real : plasma effective charge + This function calculates the bootstrap current fraction, + using the Nevins et al method, 4/11/90. + AEA FUS 251: A User's Guide to the PROCESS Systems Code + """ + # Calculate peak electron beta + + betae0 = ( + physics_variables.ne0 + * physics_variables.te0 + * 1.0e3 + * constants.echarge + / (bt**2 / (2.0 * constants.rmu0)) + ) + + # Call integration routine - if physics_variables.ips == 0: - po.ocmmnt( - self.outfile, " Pfirsch-Schluter current fraction not calculated" - ) - elif physics_variables.ips == 1: - po.ocmmnt( - self.outfile, - " (SCENE Pfirsch-Schluter current fraction scaling used)", - ) + ainteg, _ = integrate.quad( + lambda y: bsinteg( + y, + dene, + ten, + bt, + rminor, + rmajor, + zeff, + alphat, + alphan, + q0, + q95, + betat, + constants.echarge, + constants.rmu0, + ), + 0, + 0.999, + epsabs=0.001, + epsrel=0.001, + ) - po.ovarrf( - self.outfile, - "Bootstrap fraction (enforced)", - "(bootipf.)", - current_drive_variables.bootipf, - "OP ", - ) - po.ovarrf( - self.outfile, - "Diamagnetic fraction (enforced)", - "(diaipf.)", - current_drive_variables.diaipf, - "OP ", - ) - po.ovarrf( - self.outfile, - "Pfirsch-Schlueter fraction (enforced)", - "(psipf.)", - current_drive_variables.psipf, - "OP ", - ) + # Calculate bootstrap current and fraction - po.ovarre( - self.outfile, - "Loop voltage during burn (V)", - "(vburn)", - physics_variables.plascur - * physics_variables.rplas - * physics_variables.facoh, - "OP ", - ) - po.ovarre( - self.outfile, - "Plasma resistance (ohm)", - "(rplas)", - physics_variables.rplas, - "OP ", - ) + aibs = 2.5 * betae0 * rmajor * bt * q95 * ainteg + return 1.0e6 * aibs / plascur - po.ovarre( - self.outfile, - "Resistive diffusion time (s)", - "(res_time)", - physics_variables.res_time, - "OP ", - ) - po.ovarre( - self.outfile, - "Plasma inductance (H)", - "(rlp)", - physics_variables.rlp, - "OP ", - ) - po.ovarrf( - self.outfile, - "Coefficient for sawtooth effects on burn V-s requirement", - "(csawth)", - physics_variables.csawth, - ) + @staticmethod + def bootstrap_fraction_sauter(): + """Bootstrap current fraction from Sauter et al scaling + author: P J Knight, CCFE, Culham Science Centre + None + This function calculates the bootstrap current fraction + using the Sauter, Angioni and Lin-Liu scaling. +

The code was supplied by Emiliano Fable, IPP Garching + (private communication). + O. Sauter, C. Angioni and Y. R. Lin-Liu, + Physics of Plasmas 6 (1999) 2834 + O. Sauter, C. Angioni and Y. R. Lin-Liu, (ERRATA) + Physics of Plasmas 9 (2002) 5140 + """ + NR = 200 - po.osubhd(self.outfile, "Fuelling :") - po.ovarre( - self.outfile, - "Ratio of He and pellet particle confinement times", - "(tauratio)", - physics_variables.tauratio, - ) - po.ovarre( - self.outfile, - "Fuelling rate (nucleus-pairs/s)", - "(qfuel)", - physics_variables.qfuel, - "OP ", - ) - po.ovarre( - self.outfile, - "Fuel burn-up rate (reactions/s)", - "(rndfuel)", - physics_variables.rndfuel, - "OP ", - ) - po.ovarrf( - self.outfile, - "Burn-up fraction", - "(burnup)", - physics_variables.burnup, - "OP ", - ) + roa = np.arange(1, NR + 1, step=1) / NR - if any(numerics.icc == 78): - po.osubhd(self.outfile, "Reinke Criterion :") - po.ovarin( - self.outfile, - "index of impurity to be iterated for divertor detachment", - "(impvardiv)", - reinke_variables.impvardiv, - ) - po.ovarre( - self.outfile, - "Minimum Impurity fraction from Reinke", - "(fzmin)", - reinke_variables.fzmin, - "OP ", + rho = np.sqrt(physics_variables.xarea / np.pi) * roa + sqeps = np.sqrt(roa * (physics_variables.rminor / physics_variables.rmajor)) + + ne = 1e-19 * np.vectorize( + lambda r: profiles_module.nprofile( + r, + physics_variables.rhopedn, + physics_variables.ne0, + physics_variables.neped, + physics_variables.nesep, + physics_variables.alphan, ) - po.ovarre( - self.outfile, - "Actual Impurity fraction", - "(fzactual)", - reinke_variables.fzactual, + )(roa) + ni = (physics_variables.dnitot / physics_variables.dene) * ne + tempe = np.vectorize( + lambda r: profiles_module.tprofile( + r, + physics_variables.rhopedt, + physics_variables.te0, + physics_variables.teped, + physics_variables.tesep, + physics_variables.alphat, + physics_variables.tbeta, ) + )(roa) + tempi = (physics_variables.ti / physics_variables.te) * tempe - def igmarcal(self): - """Routine to calculate ignition margin - author: P J Knight, CCFE, Culham Science Centre - outfile : input integer : Fortran output unit identifier - This routine calculates the ignition margin at the final point - with different scalings. - AEA FUS 251: A User's Guide to the PROCESS Systems Code - """ + zef = np.full_like(tempi, physics_variables.zeff) # Flat Zeff profile assumed - po.oheadr(self.outfile, "Energy confinement times, and required H-factors :") - po.ocmmnt( - self.outfile, - f"{'':>5}{'scaling law':<25}{'confinement time (s)':<25}H-factor for", - ) - po.ocmmnt( - self.outfile, - f"{'':>34}{'for H = 1':<20}power balance", + # mu = 1/safety factor + # Parabolic q profile assumed + mu = 1 / ( + physics_variables.q0 + + (physics_variables.q - physics_variables.q0) * roa**2 ) + amain = np.full_like(mu, physics_variables.afuel) + zmain = np.full_like(mu, 1.0 + physics_variables.fhe3) - for iisc in range(32, 48): - ( - physics_variables.kappaa, - ptrez, - ptriz, - taueez, - taueiz, - taueffz, - powerhtz, - ) = self.pcond( - physics_variables.afuel, - physics_variables.palpmw, - physics_variables.aspect, - physics_variables.bt, - physics_variables.dnitot, - physics_variables.dene, - physics_variables.dnla, - physics_variables.eps, - 1.0, - physics_variables.iinvqd, - iisc, - physics_variables.ignite, - physics_variables.kappa, - physics_variables.kappa95, - physics_variables.pchargemw, - current_drive_variables.pinjmw, - physics_variables.plascur, - physics_variables.pcoreradpv, - physics_variables.rmajor, - physics_variables.rminor, - physics_variables.te, - physics_variables.ten, - physics_variables.tin, - physics_variables.q, - physics_variables.qstar, - physics_variables.vol, - physics_variables.xarea, - physics_variables.zeff, - ) + if ne[NR - 1] == 0.0: + ne[NR - 1] = 1e-4 * ne[NR - 2] + ni[NR - 1] = 1e-4 * ni[NR - 2] - physics_variables.hfac[iisc - 1] = self.fhfac(iisc) + if tempe[NR - 1] == 0.0: + tempe[NR - 1] = 1e-4 * tempe[NR - 2] + tempi[NR - 1] = 1e-4 * tempi[NR - 2] - po.ocmmnt( - self.outfile, - f"{'':>2}{f2py_compatible_to_string(physics_variables.tauscl[iisc-1]):<32}" - f"{taueez:<26.3f}{physics_variables.hfac[iisc - 1]:.3f}", + # Calculate total bootstrap current (MA) by summing along profiles + # Looping from 2 because dcsa etc should return 0 @ j == 1 + nr_rng = np.arange(2, NR) + nr_rng_1 = nr_rng - 1 + + drho = rho[nr_rng] - rho[nr_rng_1] + da = 2 * np.pi * rho[nr_rng_1] * drho # area of annulus + dlogte_drho = (np.log(tempe[nr_rng]) - np.log(tempe[nr_rng_1])) / drho + dlogti_drho = (np.log(tempi[nr_rng]) - np.log(tempi[nr_rng_1])) / drho + dlogne_drho = (np.log(ne[nr_rng]) - np.log(ne[nr_rng_1])) / drho + + jboot = ( + 0.5 + * ( + dcsa( + nr_rng, + NR, + physics_variables.rmajor, + physics_variables.bt, + physics_variables.triang, + ne, + ni, + tempe, + tempi, + mu, + rho, + zef, + sqeps, + ) + * dlogne_drho + + hcsa( + nr_rng, + NR, + physics_variables.rmajor, + physics_variables.bt, + physics_variables.triang, + ne, + ni, + tempe, + tempi, + mu, + rho, + zef, + sqeps, + ) + * dlogte_drho + + xcsa( + nr_rng, + NR, + physics_variables.rmajor, + physics_variables.bt, + physics_variables.triang, + mu, + sqeps, + tempi, + tempe, + amain, + zmain, + ni, + ne, + rho, + zef, + ) + * dlogti_drho + ) + * 1.0e6 + * ( + -physics_variables.bt + / (0.2 * np.pi * physics_variables.rmajor) + * rho[nr_rng_1] + * mu[nr_rng_1] ) + ) # A/m2 + + return np.sum(da * jboot, axis=0) / physics_variables.plascur def fhfac(self, is_): """Function to find H-factor for power balance @@ -4729,6 +4694,7 @@ def fhfac(self, is_): using the given energy confinement scaling law. AEA FUS 251: A User's Guide to the PROCESS Systems Code """ + physics_module.iscz = is_ return root_scalar(self.fhz, bracket=(0.01, 150), xtol=0.003).root @@ -4743,6 +4709,7 @@ def fhz(self, hhh): value of hhh, the confinement time H-factor. AEA FUS 251: A User's Guide to the PROCESS Systems Code """ + ( physics_variables.kappaa, ptrez, @@ -4796,19 +4763,19 @@ def fhz(self, hhh): # calculation (i.e. whether device is ignited) if physics_variables.ignite == 0: - fhz = fhz - current_drive_variables.pinjmw / physics_variables.vol + fhz -= current_drive_variables.pinjmw / physics_variables.vol # Include the radiation power if requested if physics_variables.iradloss == 0: - fhz = fhz + physics_variables.pradpv + fhz += physics_variables.pradpv elif physics_variables.iradloss == 1: - fhz = fhz + physics_variables.pcoreradpv + fhz += physics_variables.pcoreradpv return fhz + @staticmethod def pcond( - self, afuel, palpmw, aspect, @@ -4893,6 +4860,7 @@ def pcond( Menard 2019, Phil. Trans. R. Soc. A 377:20170440 Kaye et al. 2006, Nucl. Fusion 46 848 """ + eps2 = eps / 2.0e0 str5 = 2.0e0 / (1.0e0 + (kappa**2)) ck2 = (0.66e0 + (1.88e0 * (np.sqrt(eps2))) - (1.54e0 * eps2)) * ( diff --git a/tests/unit/test_physics.py b/tests/unit/test_physics.py index 18cf72d0d6..dfc8b50251 100644 --- a/tests/unit/test_physics.py +++ b/tests/unit/test_physics.py @@ -2,6 +2,7 @@ from typing import Any, NamedTuple from process.fortran import ( + constants, physics_variables, physics_module, current_drive_variables, @@ -10,7 +11,18 @@ ) import numpy import pytest -from process.physics import Physics +from process.physics import ( + Physics, + bpol, + diamagnetic_fraction_scene, + diamagnetic_fraction_hender, + ps_fraction_scene, + plasc, + culblm, + conhas, + vscalc, + rether, +) from process.plasma_profiles import PlasmaProfile from process.current_drive import CurrentDrive from process.impurity_radiation import initialise_imprad @@ -26,26 +38,26 @@ def physics(): return Physics(PlasmaProfile(), CurrentDrive()) -def test_diamagnetic_fraction_hender(physics): +def test_diamagnetic_fraction_hender(): """Test diamagnetic_fraction_hender().""" beta = 0.14 - diacf = physics.diamagnetic_fraction_hender(beta) + diacf = diamagnetic_fraction_hender(beta) assert diacf == pytest.approx(0.05, abs=0.0001) -def test_diamagnetic_fraction_scene(physics): +def test_diamagnetic_fraction_scene(): """Test diamagnetic_fraction_scene.""" beta = 0.15 q95 = 3.0 q0 = 1.0 - diacf = physics.diamagnetic_fraction_scene(beta, q95, q0) + diacf = diamagnetic_fraction_scene(beta, q95, q0) assert diacf == pytest.approx(0.0460, abs=0.0001) -def test_ps_fraction_scene(physics): +def test_ps_fraction_scene(): """Test ps_fraction_scene.""" beta = 0.15 - pscf = physics.ps_fraction_scene(beta) + pscf = ps_fraction_scene(beta) assert pscf == pytest.approx(-0.0135, abs=0.0001) @@ -650,8 +662,8 @@ def test_culcur(culcurparam, monkeypatch, physics): ), ), ) -def test_plasc(arguments, expected, physics): - assert physics.plasc(**arguments) == pytest.approx(expected) +def test_plasc(arguments, expected): + assert plasc(**arguments) == pytest.approx(expected) @pytest.mark.parametrize( @@ -667,6 +679,7 @@ def test_plasc(arguments, expected, physics): "kappa": 1.85, "delta": 0.5, "perim": 24, + "rmu0": constants.rmu0, }, 3.4726549397470703, ), @@ -680,6 +693,7 @@ def test_plasc(arguments, expected, physics): "kappa": 1.85, "delta": 0.5, "perim": 24, + "rmu0": constants.rmu0, }, 2.958739919272374, ), @@ -693,21 +707,22 @@ def test_plasc(arguments, expected, physics): "kappa": 1.85, "delta": 0.5, "perim": 24, + "rmu0": constants.rmu0, }, 0.8377580413333333, ), ), ) -def test_bpol(arguments, expected, physics): - assert physics.bpol(**arguments) == pytest.approx(expected) +def test_bpol(arguments, expected): + assert bpol(**arguments) == pytest.approx(expected) -def test_culblm(physics): - assert physics.culblm(12, 4.879, 18300000, 2.5) == pytest.approx(0.0297619) +def test_culblm(): + assert culblm(12, 4.879, 18300000, 2.5) == pytest.approx(0.0297619) -def test_conhas(physics): - assert physics.conhas(5, 5, 12, 0.5, 0.33, 1.85, 2e3) == pytest.approx( +def test_conhas(): + assert conhas(5, 5, 12, 0.5, 0.33, 1.85, 2e3, constants.rmu0) == pytest.approx( 2.518876726889116 ) @@ -1349,7 +1364,7 @@ class VscalcParam(NamedTuple): ), ), ) -def test_vscalc(vscalcparam, physics): +def test_vscalc(vscalcparam): """ Automatically generated Regression Unit Test for vscalc. @@ -1359,7 +1374,7 @@ def test_vscalc(vscalcparam, physics): :type vscalcparam: vscalcparam """ - phiint, rlp, vsbrn, vsind, vsres, vsstt = physics.vscalc( + phiint, rlp, vsbrn, vsind, vsres, vsstt = vscalc( csawth=vscalcparam.csawth, eps=vscalcparam.eps, facoh=vscalcparam.facoh, @@ -1371,6 +1386,7 @@ def test_vscalc(vscalcparam, physics): rplas=vscalcparam.rplas, tburn=vscalcparam.tburn, t_fusion_ramp=vscalcparam.t_fusion_ramp, + rmu0=constants.rmu0, ) assert phiint == pytest.approx(vscalcparam.expected_phiint) @@ -1518,8 +1534,8 @@ def test_phyaux(phyauxparam, monkeypatch, physics): assert taup == pytest.approx(phyauxparam.expected_taup) -def test_rether(physics): - assert physics.rether( +def test_rether(): + assert rether( 1.0, 1.45, 7.5e19, 17.81065204, 12.0, 13.0, 0.43258985 ) == pytest.approx(0.028360489673597476)