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{split}\begin{align}
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{align}\end{split}\]
The result is a curious system which depends on both the temperatures and the
temperature gradients:
\[\begin{split}\begin{align}
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{align}\end{split}\]
where
\[\begin{equation}
z_i \doteq - \frac{a}{T_i} \frac{\partial T_i}{\partial r}
\quad\mbox{and}\quad
z_e \doteq - \frac{a}{T_e} \frac{\partial T_e}{\partial r}
\end{equation}\]
It is important to note the connection between profiles and gradients. Specifically, if we enforce the following pedestal
boundary conditions at \(r=r_*\):
\[\begin{equation}
T_\sigma(r_*) = T_\sigma^* \; .
\end{equation}\]
Then the gradients \(z_\sigma\) uniquely determine the temperature profiles, \(T_\sigma\):
\[\begin{equation}
T_\sigma(r) = T_\sigma^*
\exp\left( \int_r^{r_*} dx \, z_\sigma(x) \right) \; .
\end{equation}\]
Formulation on a discrete grid
On a discrete grid \(r_j\), the temperature profile can be approximately determined using
the trapezoidal rule
\[\begin{equation}
T_\sigma(r_{j-1}) = T_\sigma(r_j) \exp \left\{
\left[ \frac{z_\sigma(r_j)+z_\sigma(r_{j-1})}{2} \right]
\left[ r_j-r_{j-1} \right] \right\} \; .
\end{equation}\]
To put the problem into discrete form, we define a vector of
independent variables (gradients) and functions (fluxes):
\[\begin{split}\begin{align}
z_{\sigma,j} = &~ z_\sigma(r_j) \; , \\
Q_{\sigma,j} = &~ Q_\sigma(r_j) \; , \\
Q^T_{\sigma,j} = &~ Q^T_\sigma(r_j) \; .
\end{align}\end{split}\]
Then, the equations to be solved are
\[\begin{equation}
\widehat{Q}_{\sigma,j} = \widehat{Q}^T_{\sigma,j} \; .
\end{equation}\]
where a hat denotes gyroBohm normalization:
\[\begin{equation}
\widehat{Q} \doteq \frac{Q}{Q_{\rm GB}}
\quad \text{where} \quad Q_{\rm GB} = n_e T_e c_s (\rho_s/a)^2 \; .
\end{equation}\]
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:
\[\widehat{Q}_{\sigma,j}(z^0) - \widehat{Q}^T_{\sigma,j}(z^0)
+ \frac{\partial \widehat{Q}_{\sigma,j}}{\partial z_{\sigma^\prime,j}} \,\delta z_{\sigma^\prime,j}
- \frac{\partial \widehat{Q}^T_{\sigma,j}}{\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^0_{\sigma,j}\}\). This can be written in terms of Jacobian matrices as
(1)\[\widehat{{\cal J}}_{{\sigma\sigma^\prime,jj^\prime}} \, \delta z_{\sigma^\prime,j^\prime} =
-\left[ \widehat{Q}_{\sigma,j}(z^0) - \widehat{Q}^T_{\sigma,j}(z^0) \right] \eta_{\sigma,j} \; ,\]
where
\[\widehat{{\cal J}}_{\sigma\sigma^\prime,jj^\prime} \doteq {\cal J}_{\sigma\sigma^\prime,jj} \delta_{jj^\prime}
-{\cal J}^T_{{\sigma\sigma^\prime,jj^\prime}} \; ,\]
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{split}\begin{equation}
{\cal J}_{{\sigma\sigma^\prime,jj^\prime}} = \begin{pmatrix}
\displaystyle \frac{\partial \widehat{Q}_{i,1}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}_{i,1}}{\partial z_{e,1}} &
0 & 0 & 0 & 0 ~ \\
\displaystyle \frac{\partial \widehat{Q}_{e,1}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}_{e,1}}{\partial z_{e,1}} &
0 & 0 & 0 & 0 ~ \\
0 & 0 &
\displaystyle \frac{\partial \widehat{Q}_{i,2}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}_{i,2}}{\partial z_{e,2}} &
0 & 0 ~ \\
0 & 0 &
\displaystyle \frac{\partial \widehat{Q}_{e,2}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}_{e,2}}{\partial z_{e,2}} &
0 & 0 ~ \\
0 & 0 & 0 & 0 &
\displaystyle \frac{\partial \widehat{Q}_{i,3}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}_{i,3}}{\partial z_{e,3}} &
~ \\
0 & 0 & 0 & 0 &
\displaystyle \frac{\partial \widehat{Q}_{e,3}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}_{e,3}}{\partial z_{e,3}} &
\end{pmatrix}
\end{equation}\end{split}\]
\[\begin{split}\begin{equation}
{\cal J}^T_{{\sigma\sigma^\prime,jj^\prime}} = \begin{pmatrix}
\displaystyle \frac{\partial \widehat{Q}^T_{i,1}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,1}}{\partial z_{e,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,1}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,1}}{\partial z_{e,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,1}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,1}}{\partial z_{e,3}}
~ \\
\displaystyle \frac{\partial \widehat{Q}^T_{e,1}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,1}}{\partial z_{e,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,1}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,1}}{\partial z_{e,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,1}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,1}}{\partial z_{e,3}}
~ \\
\displaystyle \frac{\partial \widehat{Q}^T_{i,2}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,2}}{\partial z_{e,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,2}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,2}}{\partial z_{e,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,2}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,2}}{\partial z_{e,3}}
~ \\
\displaystyle \frac{\partial \widehat{Q}^T_{e,2}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,2}}{\partial z_{e,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,2}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,2}}{\partial z_{e,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,2}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,2}}{\partial z_{e,3}}
~ \\
\displaystyle \frac{\partial \widehat{Q}^T_{i,3}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,3}}{\partial z_{e,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,3}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,3}}{\partial z_{e,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,3}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}^T_{i,3}}{\partial z_{e,3}}
~ \\
\displaystyle \frac{\partial \widehat{Q}^T_{e,3}}{\partial z_{i,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,3}}{\partial z_{e,1}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,3}}{\partial z_{i,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,3}}{\partial z_{e,2}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,3}}{\partial z_{i,3}} &
\displaystyle \frac{\partial \widehat{Q}^T_{e,3}}{\partial z_{e,3}}
\end{pmatrix}
\end{equation}\end{split}\]
An important quantity to measure after a Newton iteration is
the residual
\[R^1_{\sigma,j} = \frac{\left[\widehat{Q}_{\sigma,j}(z^1)-\widehat{Q}^T_{\sigma,j}(z^1)\right]^2}{
\left[\widehat{Q}_{\sigma,j}(z^1)\right]^2+\left[\widehat{Q}^T_{\sigma,j}(z^1)\right]^2}\]
If, after a Newton step, any \(R^1_{\sigma,j} > R^0_{\sigma,j}\) 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 \(\widehat{Q}\) was included, there would be no fixed-point
iteration component.