7. The BESSER formal solver
To solve the radiation transfer equations, PORTA applies the generalization of the short-characteristics (SC) method of Kunasz, P. and Auer, L.H. (1988) [1] to the polarized case. For a given frequency and direction, we consider three consecutive spatial points \(M\), \(O\), and \(P\) along the propagation direction, with \(M\) the upwind points, \(P\) the downwind points, and \(O\) the point where the Stokes parameters need to be computed.
Fig. 7.1 Short-characteristics in a three-dimensional Cartesian rectilinear grid.
The aim of the original SC method is to solve the integral form of the radiation transfer equation:
where \(I_O\) is the intensity we want to compute, \(I_M\) the intensity in the upwind point, \(S\left(t\right)\) is the source function at \(\tau = t\), and \(\tau_{MO}\) is the optical distance between points \(M\) and \(O\). In our case, the equation to solve is:
where \(\vec{S}_O\) is the vector of Stokes parameters we want to compute, \(\vec{S}_M\) is the Stokes vector in the upwind points, and \(S_{\rm eff}\left(t\right)\) is the effective source function given by:
with \(\vec{K}'\left(t\right) = \vec{K}\left(t\right)/\eta_I - \vec{1}.\)
The original SC method is based on the approximation of parabolic interpolation of the source function between the upwind point \(M\), the grid point \(O\), and the downwind point \(P\) (see Fig. 7.1). In 2D and 3D grids, the upwind and downwind points of the SC do not generally coincide with any spatial grid node and the radiation transfer quantities (i.e., the emission and absorption coefficients) have to be interpolated from the nearby 9-point (biquadratic case) or 4-point (bilinear case) stencils of the discrete grid points. Bilinear interpolation is sufficient in the fine grids of today’s MHD models and this option is implemented in PORTA. The biquadratic interpolation can lead to more accurate solution but it is significantly more time consuming. The upwind Stokes vector \(\vec{S}_M\) needs to be interpolated from the same grid nodes. Proper topological ordering of the grid points is therefore necessary for every direction of the short characteristics. The intersection points \(M\) and \(P\) may be located on a vertical plane of the grid instead of a horizontal one.
To avoid the overshooting problems of the parabolic interpolation, Štěpán, J., and Trujillo Bueno, J. (2013) [2] developed an accurate formal solver, dubbed BESSER (BEzier Spline SolvER). The algorithm implements an interpolation scheme based on the use of monotonic Bézier splines (see Auer, L. (2003) [3] for an alternative implementation with splines), imposing a continuous first derivative of the interpolant.
The control points of the splines are used to preserve monotonicity of the interpolant (it is contained in an envelope defined by the tangents of the spline in its endpoints and of the control point which is located in the intersection of these tangents, Fig. 7.2). By imposing a continuous first derivative of the interpolant, we achieve a smooth connection of the Bézier splines in the central point \(O\). This implementation leads to a symmetrical interpolation independently of the choice of the interpolation direction (\(MOP\) for one direction of the ray propagation or \(POM\) for the opposite direction of the ray). In addition, our BESSER method always provides reliable values for the diagonal of the \(\Lambda\)-operator, i.e., in the interval [0,1), used in methods based on the Jacobi iteration. See Fig. 7.2 for a comparison between the parabolic and BESSER implementations.
Fig. 7.2 Parabolic and BESSER interpolation using the three successive points M, O, and P. Dotted line: parabolic interpolation. Solid line: interpolation using our BESSER method with continuous derivative at point O.
Given a quantity \(y\) (e.g., the source function) defined at three successive points \(x_M\), \(x_O\), andi \(x_P\), we use two quadratic Bézier splines to interpolate \(y\) between points \(M\) and \(O\) and between points \(O\) and \(P\) (see Fig. 7.2). First, we look for an optimal interpolation in the interval \(MO\). For the sake of simplicity, we parametrize the \(x\)-coordinate in this interval by a dimensionless parameter \(u=(x-x_M)/h_M\), where \(h_M=x_O-x_M\). The Bézier spline in the interval \(MO\) is a parabola passing through points \(M\) and \(O\). The derivatives at such points are defined by the position of the control point whose \(y\)-coordinate is \(c_M\) (see Fig. 7.2). The equation for such a spline reads (Auer, L. 2003 [3]):
Similarly, one can define a Bézier spline between points \(O\) and \(P\) by doing the formal changes \(y_M\rightarrow y_O\), \(y_O\rightarrow y_P\), and \(u=(x-x_O)/h_P\), where \(h_P=x_P-x_O\); the \(y\)-coordinate of the ensuing control point is denoted by \(c_P\) (see Fig. 7.2). We look for the values of \(c_M\) and \(c_P\) that satisfy the following conditions: (1) if the sequence \(y_M\), \(y_O\), \(y_P\) is monotonic, then the interpolation is monotonic in the whole interval \([x_M,x_P]\); (2) if the sequence of \(y_i\) values is not monotonic, then the interpolant has the only local extremum at \(O\); and (3) the first derivative of the interpolant at point \(O\) should be continuous. The ensuing algorithm proceeds as follows:
Calculate the quantities \(d_M=(y_O-y_M)/h_M\), \(d_P=(y_P-y_O)/h_P\).
If the sequence \(y_M\), \(y_O\), \(y_P\) is not monotonic (i.e., if \(d_Md_P\le 0\)), then set \(c_M=c_P=y_O\) and exit the algorithm. The derivative of the splines at point \(O\) is equal to zero, leading to a local extremum at the central point.
Estimate the derivative at point \(O\),
\({y'}_O = \frac{h_Md_P + h_Pd_M}{h_m+h_p}\)
(see Eq. (7) of Auer, L. (2003) [3], and references therein). This derivative is equal to that provided at the same point by parabolic interpolation among points \(M\), \(O\), and \(P\). Moreover, in contrast with Eq. (12) of Auer, L. (2003) [3], it is an expression that relates \(y'\) (e.g., the source function derivative) linearly with the \(y\)-values (e.g., with the source function values).
Calculate the initial positions of the control points, \(c_M=y_O-\frac{h_M}{2}{y'}_O\), and \(c_P=y_O+\frac{h_P}{2}{y'}_O\). The control points calculated this way lead to a unique parabolic interpolation among points \(MOP\). If the algorithm is stopped here, the resulting formal solver would be equivalent to the standard parabolic interpolation.
Check that \(min(y_M,y_O)\le c_M\le max(y_M,y_O)\). If the condition is satisfied, then go to step (7), otherwise continue with step (6).
If the condition in step (5) is not satisfied, then there is an overshoot of the interpolant in the interval \(MO\). Set \(c_M=y_M\), so that the first derivative at \(M\) is equal to zero and the overshoot is corrected. Since the value of \(c_P\) is not of interest for the formal solution between \(M\) and \(O\), exit the algorithm.
Check if \(min(y_O,y_P)\le c_P\le max(y_O,y_P)\). If this conditionis not satisfied, then set \(c_P=y_P\) so that the overshoot in the interval \(OP\) is corrected.
Calculate the new derivative at \(O\), \({y'}_O=(c_P-y_O)/(h_P/2)\), using the corrected value of \(c_P\) calculated in step (7).
Calculate a new \(c_M\) value to keep the derivative at \(O\) smooth. It is easy to realize that this change cannot produce an overshoot in the \(MO\) interval, hence the solution remains monotonic with a continuous derivative.
Steps (8) and (9) of the above-mentioned algorithm deal with the correction of the overshoots in the downwind interval followed by modification of the \(c_M\) upwind control point value. We have found that it is suitable to guarantee the smoothness of the derivative, and this can be done with only a small increase in the computing time with respect to the parabolic method. Our BESSER interpolation is stable; that is, the interpolant varies smoothly with smooth changes of the \(M\), \(O\), and \(P\) points. No abrupt changes of the splines occur that could negatively affect the stability of the iterative method.
In contrast to some other formal solvers based on the idea of quadratic Bézier splines, our BESSER algorithm guarantees that a monotonic sequence of the \(MOP\) points leads to a monotonic interpolant in all situations. This fact is of critical importance in 2D and 3D grids in which \(\tau_{MO}\) and \(\tau_{OP}\) may differ significantly because of unpredictable intersections of the grid planes, especially if periodic boundary conditions are considered. Such large differences often lead to overshoots, unless treated by BESSER or a similarly suitable strategy.