TGYRO solver algorithm

TGYRO solves the steady-state transport problem; that is, the transport equations with ${\partial }_t\to 0$. This algorithm was originally envisioned to form the basis for a time-implicit solver, but frankly there has been little interest or need for this capability. Instead, flexibility in dealing with multiple ions and solver robustness have been emphasized.

Formulation

To describe the algorithm, we will restrict attention to coupled $T_i$ and $T_e$ evolution, even though density and rotation can also be evolved.

$$\begin{array}{r}\begin{array}{rl}\frac{1}{V^{\prime }(r)}\frac{\partial }{\partial r}[V^{\prime }(r)Q_i(r)]=& S_i,\\ \frac{1}{V^{\prime }(r)}\frac{\partial }{\partial r}[V^{\prime }(r)Q_e(r)]=& S_e.\end{array}\end{array}$$

Here, the energy fluxes are the taken to be the sum of neoclassical and turbulent transport:

$$\begin{array}{r}\begin{array}{rl}{Q}_i=& Q_i^{\mathrm{N}\mathrm{e}\mathrm{o}}+Q_i^{\mathrm{T}\mathrm{u}\mathrm{r}\mathrm{b}}\\ Q_e=& Q_e^{\mathrm{N}\mathrm{e}\mathrm{o}}+Q_e^{\mathrm{T}\mathrm{u}\mathrm{r}\mathrm{b}}\end{array}\end{array}$$

The total ion and electron sources, $S_e$ and $S_i$, are described in more detail in the Equations solved in TGYRO .

Some comments regarding units

In TGYRO, we have found it convenient to use CGS units rather than employing some variant of the more popular dimensionless normalizations. Thus, we have

$$\begin{array}{r}\begin{array}{rl}\mathrm{S}\mathrm{o}\mathrm{u}\mathrm{r}\mathrm{c}\mathrm{e}:& S\sim \frac{\mathrm{e}\mathrm{r}\mathrm{g}}{{\mathrm{c}\mathrm{m}}^3\mathrm{s}}\\ \mathrm{E}\mathrm{n}\mathrm{e}\mathrm{r}\mathrm{g}\mathrm{y}\mathrm{F}\mathrm{l}\mathrm{u}\mathrm{x}:& Q\sim \frac{\mathrm{e}\mathrm{r}\mathrm{g}}{{\mathrm{c}\mathrm{m}}^2\mathrm{s}}\\ \mathrm{P}\mathrm{o}\mathrm{w}\mathrm{e}\mathrm{r}:& P\sim \frac{\mathrm{e}\mathrm{r}\mathrm{g}}{\mathrm{s}}\to {\int }_0^{r}dx{V}^{\prime }(x)S(x)\end{array}\end{array}$$

Solution strategy

Rather than solving the equations directly, we prefer to solve the volume-integrated form of the equation so that we can deal directly with the fluxes:

$$\begin{array}{r}\begin{array}{rl}{Q}_i^T(r)\doteq \frac{1}{V^{\prime }(r)}& {\int }_0^{r}dx{V}^{\prime }(x)S_i\\ Q_e^T(r)\doteq \frac{1}{V^{\prime }(r)}& {\int }_0^{r}dx{V}^{\prime }(x)S_e\end{array}\end{array}$$

The result is a curious system which depends on both the temperatures and the temperature gradients:

$$\begin{array}{r}\begin{array}{rl}{Q}_i(z_i,z_e,T_i,T_e)-Q_i^T(T_i,T_e)=& 0\\ Q_e(z_i,z_e,T_i,T_e)-Q_e^T(T_i,T_e)=& 0\end{array}\end{array}$$

where

$$z_i\doteq -\frac{a}{T_i}\frac{\partial T_i}{\partial r}\text{and}{z}_e\doteq -\frac{a}{T_e}\frac{\partial T_e}{\partial r}$$

It is important to note the connection between profiles and gradients. Specifically, if we enforce the following pedestal boundary conditions at $r=r_{\ast }$:

$$T_{\sigma }(r_{\ast })=T_{\sigma }^{\ast }.$$

Then the gradients $z_{\sigma }$ uniquely determine the temperature profiles, $T_{\sigma }$:

$$T_{\sigma }(r)=T_{\sigma }^{\ast }\exp ({\int }_r^{r_{\ast }}dx{z}_{\sigma }(x)).$$

Formulation on a discrete grid

On a discrete grid $r_j$, the temperature profile can be approximately determined using the trapezoidal rule

$$T_{\sigma }(r_{j-1})=T_{\sigma }(r_j)\exp \left\{\left[\frac{z_{\sigma }(r_j)+z_{\sigma }(r_{j-1})}{2}\right][r_j-r_{j-1}]\right\}.$$

To put the problem into discrete form, we define a vector of independent variables (gradients) and functions (fluxes):

$$\begin{array}{r}\begin{array}{rl}{z}_{\sigma ,j}=& z_{\sigma }(r_j),\\ Q_{\sigma ,j}=& Q_{\sigma }(r_j),\\ Q_{\sigma ,j}^T=& Q_{\sigma }^T(r_j).\end{array}\end{array}$$

Then, the equations to be solved are

$${\hat{Q}}_{\sigma ,j}={\hat{Q}}_{\sigma ,j}^T.$$

where a hat denotes gyroBohm normalization:

$$\hat{Q}\doteq \frac{Q}{Q_{\mathrm{G}\mathrm{B}}}\text{where}{Q}_{\mathrm{G}\mathrm{B}}=n_e{T}_e{c}_s({\rho }_s/a{)}^2.$$

The goal is to apply Newton’s method in a way which is as accurate as possible while still minimizing evaluation of the expensive functions $Q_{\sigma ,j}$. Operationally, we make the key assumption that the transport fluxes depend only locally on the gradients (which is approximately true when quantities are normalzied to the gyroBohm unit of flux), so that the Jacobian associated with $Q_{\sigma ,j}$ is block diagonal:

$${\hat{Q}}_{\sigma ,j}(z^0)-{\hat{Q}}_{\sigma ,j}^T(z^0)+\frac{\partial {\hat{Q}}_{\sigma ,j}}{\partial z_{{\sigma }^{\prime },j}}\delta z_{{\sigma }^{\prime },j}-\frac{\partial {\hat{Q}}_{\sigma ,j}^T}{\partial z_{{\sigma }^{\prime },j^{\prime }}}\delta z_{{\sigma }^{\prime },j^{\prime }}=0.$$

Above, we have used the shorthand $z\doteq \{z_{\sigma ,j}\}$ and $z^0\doteq \{z_{\sigma ,j}^0\}$. This can be written in terms of Jacobian matrices as

(1)$${\hat{\mathcal{J}}}_{\sigma {\sigma }^{\prime },j{j}^{\prime }}\delta z_{{\sigma }^{\prime },j^{\prime }}=-[{\hat{Q}}_{\sigma ,j}(z^0)-{\hat{Q}}_{\sigma ,j}^T(z^0)]{\eta }_{\sigma ,j},$$

where

$${\hat{\mathcal{J}}}_{\sigma {\sigma }^{\prime },j{j}^{\prime }}\doteq {\mathcal{J}}_{\sigma {\sigma }^{\prime },jj}{\delta }_{j{j}^{\prime }}-{\mathcal{J}}_{\sigma {\sigma }^{\prime },j{j}^{\prime }}^T,$$

and the quantity $z^1=z^0+\delta z$ is the Newton update for the vector $z$. In Eq. (1), we have introduced a relaxation parameter ${\eta }_{\sigma ,j}$. Note that this method generalizes to an arbitrary number of gradients and fluxes per gridpoint. In the case of three radial gridpoints, $\{r_1,r_2,r_3\}$, the Jacobian matrices have the explicit forms

$$\begin{array}{r}{\mathcal{J}}_{\sigma {\sigma }^{\prime },j{j}^{\prime }}=\left(\begin{array}{cccccc}{\displaystyle \frac{\partial {\hat{Q}}_{i,1}}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,1}}{\partial z_{e,1}}}& 0& 0& 0& 0\\ {\displaystyle \frac{\partial {\hat{Q}}_{e,1}}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,1}}{\partial z_{e,1}}}& 0& 0& 0& 0\\ 0& 0& {\displaystyle \frac{\partial {\hat{Q}}_{i,2}}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,2}}{\partial z_{e,2}}}& 0& 0\\ 0& 0& {\displaystyle \frac{\partial {\hat{Q}}_{e,2}}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,2}}{\partial z_{e,2}}}& 0& 0\\ 0& 0& 0& 0& {\displaystyle \frac{\partial {\hat{Q}}_{i,3}}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,3}}{\partial z_{e,3}}}& \\ 0& 0& 0& 0& {\displaystyle \frac{\partial {\hat{Q}}_{e,3}}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,3}}{\partial z_{e,3}}}& \end{array}\right)\end{array}$$

$$\begin{array}{r}{\mathcal{J}}_{\sigma {\sigma }^{\prime },j{j}^{\prime }}^T=\left(\begin{array}{cccccc}{\displaystyle \frac{\partial {\hat{Q}}_{i,1}^T}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,1}^T}{\partial z_{e,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,1}^T}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,1}^T}{\partial z_{e,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,1}^T}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,1}^T}{\partial z_{e,3}}}\\ {\displaystyle \frac{\partial {\hat{Q}}_{e,1}^T}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,1}^T}{\partial z_{e,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,1}^T}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,1}^T}{\partial z_{e,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,1}^T}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,1}^T}{\partial z_{e,3}}}\\ {\displaystyle \frac{\partial {\hat{Q}}_{i,2}^T}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,2}^T}{\partial z_{e,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,2}^T}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,2}^T}{\partial z_{e,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,2}^T}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,2}^T}{\partial z_{e,3}}}\\ {\displaystyle \frac{\partial {\hat{Q}}_{e,2}^T}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,2}^T}{\partial z_{e,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,2}^T}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,2}^T}{\partial z_{e,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,2}^T}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,2}^T}{\partial z_{e,3}}}\\ {\displaystyle \frac{\partial {\hat{Q}}_{i,3}^T}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,3}^T}{\partial z_{e,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,3}^T}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,3}^T}{\partial z_{e,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,3}^T}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{i,3}^T}{\partial z_{e,3}}}\\ {\displaystyle \frac{\partial {\hat{Q}}_{e,3}^T}{\partial z_{i,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,3}^T}{\partial z_{e,1}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,3}^T}{\partial z_{i,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,3}^T}{\partial z_{e,2}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,3}^T}{\partial z_{i,3}}}& {\displaystyle \frac{\partial {\hat{Q}}_{e,3}^T}{\partial z_{e,3}}}\end{array}\right)\end{array}$$

An important quantity to measure after a Newton iteration is the residual

$$R_{\sigma ,j}^1=\frac{{[{\hat{Q}}_{\sigma ,j}(z^1)-{\hat{Q}}_{\sigma ,j}^T(z^1)]}^2}{{[{\hat{Q}}_{\sigma ,j}(z^1)]}^2+{[{\hat{Q}}_{\sigma ,j}^T(z^1)]}^2}$$

If, after a Newton step, any $R_{\sigma ,j}^1>R_{\sigma ,j}^0$ is not reduced, some strategy must be adopted to modify the gradient vector $z^1$ and/or the target. Note that there are two distinct iterations:

• A Newton iteration, which is rapidly convergent given that one is close to a root and the $qhat$ are smooth functions,

• A fixed-point iteration following the Newton iteration, because the weak profile variation of $qhat$ was ignored

If the temperature dependence of $\hat{Q}$ was included, there would be no fixed-point iteration component.

Computation of the Jacobian

We approximate the derivatives in the Jacobian matrix using a forward difference approximation

$$\frac{\partial {\hat{Q}}_{\sigma ,j}}{\partial z_{{\sigma }^{\prime },j^{\prime }}}\simeq \frac{{\hat{Q}}_{\sigma ,j}(z_{{\sigma }^{\prime },j^{\prime }}+\mathrm{\Delta }z)-{\hat{Q}}_{\sigma ,j}(z_{{\sigma }^{\prime },j^{\prime }})}{\mathrm{\Delta }z}$$

A desireable feature of this approximation is that the iteration scheme, Eq.~(ref{eq.newton}) if it converges, will converge to the exact root of the original equations without any influence of the finite-difference truncation error.