Issue 
A&A
Volume 624, April 2019



Article Number  A104  
Number of page(s)  15  
Section  Numerical methods and codes  
DOI  https://doi.org/10.1051/00046361/201834761  
Published online  18 April 2019 
A novel fourthorder WENO interpolation technique
A possible new tool designed for radiative transfer
^{1}
Istituto Ricerche Solari Locarno (IRSOL), 6605 LocarnoMonti, Switzerland
email: gioele.janett@irsol.ch
^{2}
Seminar for Applied Mathematics (SAM) ETHZ, DMATH, 8093 Zurich, Switzerland
^{3}
KiepenheuerInstitut für Sonnenphysik (KIS), Schöneckstrasse 6, 79104 Freiburg i. Br., Germany
Received:
3
December
2018
Accepted:
15
March
2019
Context. Several numerical problems require the interpolation of discrete data that present at the same time (i) complex smooth structures and (ii) various types of discontinuities. The radiative transfer in solar and stellar atmospheres is a typical example of such a problem. This calls for highorder wellbehaved techniques that are able to interpolate both smooth and discontinuous data.
Aims. This article expands on different nonlinear interpolation techniques capable of guaranteeing highorder accuracy and handling discontinuities in an accurate and nonoscillatory fashion. The final aim is to propose new techniques which could be suitable for applications in the context of numerical radiative transfer.
Methods. We have proposed and tested two different techniques. Essentially nonoscillatory (ENO) techniques generate several candidate interpolations based on different substencils. The smoothest candidate interpolation is determined from a measure for the local smoothness, thereby enabling the essentially nonoscillatory property. Weighted ENO (WENO) techniques use a convex combination of all candidate substencils to obtain highorder accuracy in smooth regions while keeping the essentially nonoscillatory property. In particular, we have outlined and tested a novel wellperforming fourthorder WENO interpolation technique for both uniform and nonuniform grids.
Results. Numerical tests prove that the fourthorder WENO interpolation guarantees fourthorder accuracy in smooth regions of the interpolated functions. In the presence of discontinuities, the fourthorder WENO interpolation enables the nonoscillatory property, avoiding oscillations. Unlike Bézier and monotonic highorder Hermite interpolations, it does not degenerate to a linear interpolation near smooth extrema of the interpolated function.
Conclusion. The novel fourthorder WENO interpolation guarantees high accuracy in smooth regions, while effectively handling discontinuities. This interpolation technique might be particularly suitable for several problems, including a number of radiative transfer applications such as multidimensional problems, multigrid methods, and formal solutions.
Key words: radiative transfer / methods: numerical
© ESO 2019
1. Introduction
In astrophysics, it is common practice to solve the radiative transfer equation in solar and stellar atmospheres by means of numerical methods. It is also known that realistic atmospheric models can be highly inhomogeneous and dynamic and often present discontinuities or sharp gradients. Interpolations are ubiquitous in such problems: for instance, they are a fundamental ingredient of multidimensional radiative transfer, of multigrid methods, and of many formal solvers. For such problems, which contain both sharp discontinuities and complex smooth features, one seeks interpolation techniques able to handle singularities in an accurate and nonoscillatory fashion and at the same time to perform with uniform highorder accuracy in smooth regions.
It is a well known fact that standard highorder interpolations tend to misrepresent the nonsmooth regions of a problem, introducing spurious oscillations near discontinuities (e.g., Richards 1991; Zhang & Martin 1997; Shu 1998). However, monotonicitypreserving strategies such as Bézier and monotonic highorder Hermite interpolations, usually sacrifice accuracy to obtain monotonicity (Aràndiga et al. 2013). In fact, near smooth extrema of the interpolated function, they degenerate to a linear interpolation, dropping to secondorder accuracy (Shu 1998). By contrast, essentially nonoscillatory (ENO) and weighted ENO (WENO) techniques accomplish the feat of approximating a function with highorder accuracy in smooth parts, while avoiding “Gibbslike” oscillations near the discontinuities (Fjordholm 2016).
Section 2 presents the ENO approximation technique and it exposes the thirdorder ENO interpolation. Section 3 presents the WENO approximation technique, with a special focus on a novel fourthorder WENO interpolation. Section 4 exposes numerical tests for different interpolations. Section 5 describes different radiative transfer problems that require highorder robust interpolation techniques. Finally, Sect. 6 provides remarks and conclusions.
2. ENO approximation
The term stencil indicates the arrangement of grid points and corresponding function values (or coefficients) that is considered for a particular interpolation or reconstruction. Fixed stencil interpolations of second or higher order of accuracy are necessarily oscillatory near a discontinuity and such oscillations do not decay in magnitude when the computational grid is refined; an explicit example will be given in Sect. 4. Such issues are well recognized and considerable efforts have been already exercised to resolve them. In a seminal paper, van Leer (1979) proposed a hybrid technique that switches from a linear polynomial (secondorder accurate) to a piecewise constant approximation (firstorder accurate) near discontinuities. Similarly, one can also make use of a quadratic polynomial interpolation and limit the polynomial degree near discontinuities. Alternatively, monotonic Hermite interpolants suppress over and undershoots by controlling the firstorder derivatives (Fritsch & Carlson 1980), whereas Bézier interpolants make use of the socalled control points to shape the curve and suppress spurious extrema, preserving monotonicity in the interpolation. These limiting strategies are very robust, produce smooth solutions, and are widely used for problems containing shocks and other discontinuities. However, such monotonicitypreserving strategies sacrifice accuracy in order to obtain monotonicity (Aràndiga et al. 2013). Their main weakness is that they necessarily degenerate to a linear interpolation near smooth extrema of the interpolated function (Shu 1998). A local extremum is any point at which the value of a function is larger (a maximum) or smaller (a minimum) than all the adjacent function values. In this context, a smooth extremum is a local extremum of the function which is at least of class C^{2}, that is, twice continuously differentiable.
By contrast, ENO techniques are highorder approximations that handle discontinuities and maintain highorder accuracy near smooth extrema. The key point of the ENO strategy is to adaptively choose the stencil in the direction of smoothness. In practice, ENO schemes generate several candidate interpolations, measure the local smoothness and choose the smoothest candidate interpolation to work with, discarding the rest. This strategy enables the essentially nonoscillatory property. There is broad and wellcited body of literature attesting to the importance and the success of this strategy: from the pioneering works by Harten et al. (1986); Harten et al. (1987) and Shu & Osher (1988); Shu & Osher (1989) to the recent survey by Zhang et al. (2016).
Most of the literature about ENO (and WENO) techniques is dedicated to reconstruction schemes and not to interpolations. However, the reconstruction of cell averages of a function is equivalent to the interpolation of the point values of its primitive. Therefore, ENO (and WENO) algorithms can be formulated as an interpolation technique (Fjordholm et al. 2011).
2.1. Divided differences
Given a set of function values {y_{k}}≔{y(x_{k})} at positions {x_{k}}, the zeroth degree divided differences are defined as
while the first degree divided differences read
In general, the jth degree forward divided differences, for j > 1, are inductively defined by
As long as the function y(x) is smooth inside [x_{k}, x_{k + j}], the standard approximation theory states that (e.g., Fjordholm 2016)
By contrast, if the function y(x) is discontinuous at some point inside [x_{k}, x_{k + j}], then^{1}
which diverges for Δx → 0. Thus, divided differences are a good measure of the smoothness of the function y(x) inside the stencil. This property plays an essential role in the adaptive choice of the stencil in ENO schemes.
2.2. ENO interpolation
Consider a function y(x) and a set of function values {y_{i}} at positions {x_{i}} with i ∈ {1, …, N}. The general pth order ENO interpolation E(x) in the interval [x_{1}, x_{N}] reads
with
and p_{i}(x) denotes the interpolating polynomial acting on the interval [x_{i}, x_{i + 1}]. This polynomial must satisfy the accuracy requirement
The choice of ENO stencil for constructing p_{i}(x) is based on the divided differences presented in Sect. 2.1. The stencil is adaptively chosen so that, starting from x_{i}, it extends in the direction where y(x) is smoothest, that is, where the divided differences are smallest.
A general ppoint stencil containing x_{i} is given by
with q ∈ {1, …, p}. For notational simplicity, the explicit dependence of the stencil on index i is not indicated. This stencil can be expanded in two ways: by adding either the left neighbor x_{i − p + q − 1} or the right neighbor x_{i + q}. One decides upon which point to add by comparing the two relevant divided differences and picking the one with a smaller absolute value. Thus, if
one adds x_{i − p + q − 1} to the stencil, otherwise, x_{i + q}.
In order to build p_{i}(x), one starts with the onepoint stencil
and adaptively adds points with the procedure described above, ending up with a ppoint stencil
We note that each substencil is a substencil of the large stencil
and always contains the considered x_{i} (but not necessarily x_{i + 1}).
Numerical analysis asserts that there is a unique polynomial q_{q}(x) of degree at most p − 1, that goes through all the data points of the substencil . Its Lagrangian expression reads
where the Lagrange basis polynomials ℓ_{k} are defined as^{2}
Hence, once the optimal candidate substencil is determined, the porder accurate polynomial p_{i}(x) in Eq. (1) is simply chosen as
with q_{q}(x) given by Eq. (4).
Even though there are very few theoretical results about the stability of ENO schemes (e.g., Fjordholm et al. 2011, 2012), these schemes are very robust and stable in practice. In fact, the ENO interpolant E(x) is essentially nonoscillatory in the sense of recovering y(x) to order 𝒪(Δx^{p}) near (but not at) the location of a discontinuity. We note that ENO schemes allow for minimal over or undershoots near discontinuities. As an explicit example of this technique, the thirdorder accurate ENO interpolation is presented in the following.
2.3. Thirdorder ENO interpolation
In order to build a thirdorder accurate p_{i}(x), one needs a threepoint stencil. Thus, starting from the onepoint stencil
one adds either the left neighbor x_{i − 1} or the right neighbor x_{i + 1}. One takes the decision by comparing the absolute values of the two relevant divided differences
The smaller one implies that the function is smoother in that stencil. Therefore, if
the twopoint stencil is taken as
otherwise,
Suppose that the twopoint stencil has been chosen. In this case, one can add either the left neighbor x_{i − 1} or the right neighbor x_{i + 2}. One takes the decision by comparing the absolute values of the two relevant divided differences
Once again, if
the threepoint stencil is taken as
otherwise,
Then, one proceeds with the Lagrange interpolation given by Eq. (4) on the selected substencil, yielding the thirdorder accurate quadratic polynomial p_{i}(x).
We note that the choice of the twopoint stencil at the previous step may lead to the threepoint stencil
In order to build p_{i}(x), the thirdorder ENO interpolation thus considers the large stencil
All the substencils , , and usually work well for globally smooth problems, while the ENO technique adaptively avoids including possible discontinuities in the selected substencil.
3. WENO approximation
WENO approximations are inspired by the ENO strategy, offering with additional advantages. Instead of using the smoothest candidate substencil to form the interpolation, WENO schemes use a convex combination of all candidate substencils to obtain highorder accuracy in smooth regions while keeping the essentially nonoscillatory property. Additionally, the WENO strategy yields smoother interpolations and avoids the use of conditional “if…else” statements, which burden the algorithm.
Liu et al. (1994) proposed this strategy, constructing the first WENO schemes. Jiang & Shu (1996) constructed arbitraryorder accurate WENO approximations for finite difference schemes and most applications use their fifthorder WENO version. General information is found in the recent surveys by Shu (2009) and Zhang et al. (2016). Many papers about WENO approximations consider their implementation to uniform grids. However, the WENO strategy extends naturally to nonuniform grids, although it becomes quite complicated, depending on the complexity of the grid structure and on the interpolation degree (ČrnjarićŽic et al. 2007).
3.1. WENO interpolation
Consider a function y(x) and a set of function values {y_{i}} at positions {x_{i}} with i ∈ {1, …, N}. The general pth order WENO interpolation W(x) in the interval [x_{1}, x_{N}] reads
where I_{i}(x) is defined in Eq. (2). Consider the ppoint stencil given by Eq. (3) to be the large stencil. The unique pthorder accurate Lagrange polynomial p(x) interpolating can be recovered as a linear combination of ℓthorder accurate Lagrange polynomials q_{m}(x) defined on a set of ℓpoint substencils with ℓ ∈ {1, …, p} and m ∈ { − p + q + ℓ, …, q}, i.e.,^{3}
where γ_{m}(x) are the socalled linear (or optimal) weights, which are always positive and univocally defined. Because of consistency one has
The key idea of the WENO strategy is to build the polynomial p_{i}(x) in Eq. (5) as a weighted combination based on Eq. (6), namely,
where the weights coefficients ω_{m}(x) must satisfy the requirements:

ω_{m}(x)≥0 ∀ m for stability;

for consistency;

ω_{m}(x)≈γ_{m}(x) ∀ m if y(x) is smooth in the large stencil ;

ω_{m}(x)≈0 if y(x) has a discontinuity in the stencil .
These (nonlinear) weights are usually defined as
where the unnormalized weights are given by
and ϵ is a small positive constant used to avoid vanishing denominators; typically ϵ = 10^{−6} (Zhang et al. 2016). A larger power a makes the weight assigned to a nonsmooth substencil approach zero faster, resulting in more dissipative WENO interpolations (Liu et al. 2018). The exponent a is a free parameter. In this work, one has either a = 1 or a = 3/2.
The nonlinear weights ω_{m}(x) rely on the smoothness indicators β_{m}, which measure the relative smoothness of the function y(x) inside the stencils . A larger β_{m} indicates a lack of smoothness of y(x) in the stencil and produces a smaller nonlinear weight ω_{m}(x). The choice of the smoothness indicators plays an essential role in the performance of WENO approximations. Within the literature, the version of Jiang & Shu (1996)
is well established. Unfortunately, indicators (7) consider uniform grids alone and are particularly dissipative (i.e., they smear sharp gradients) for loworder WENO interpolations. In the following, a thirdorder and a novel fourthorder accurate WENO interpolations are presented.
3.2. Thirdorder WENO interpolation
Compared with higherorder versions, the thirdorder WENO interpolation is more robust for the treatment of discontinuous problems, it uses fewer grid points^{4}, and it provides a suitable compromise between computational cost and accuracy (Liu et al. 2018).
The unique quadratic Lagrange polynomial p(x) interpolating the threepoint large stencil
can be written as
in other words, as a linear combination of the linear Lagrange interpolations
based, respectively, on the two twopoint substencils
The linear weights are given by
and satisfy γ_{1}(x)+γ_{2}(x)=1. The linear weights depend just on the grid geometry and not on the function values. The thirdorder WENO interpolation in the interval [x_{i}, x_{i + 1}] reads
where the nonlinear weights are defined as
with
and ϵ = 10^{−6}.
3.2.1. Smoothness indicators for uniform grids
In uniform grids, the conventional choice for the smoothness indicators (7) yields
Unfortunately, this version renders the thirdorder WENO interpolation too dissipative.
Alternatively, the mapping function of the WENOM scheme (Henrick et al. 2005) and the global stencil indicator of the WENOZ scheme (Borges et al. 2008) effectively improve the accuracy of highorder WENO schemes, but they are not satisfactory for the thirdorder one. More recently, Wu et al. (2015) proposed the global smoothness indicators for the less dissipative thirdorder WENON3 and WENONP3 schemes. Similarly, Gande et al. (2017) proposed the global smoothness indicator of the WENOF3 scheme. However, such methods are usually not satisfactory in terms of accuracy or generate apparent oscillatory interpolations.
For these reasons, Liu et al. (2018) proposed smoothness indicators that make use of all the three points of the stencil, namely,
Numerical tests demonstrate that these new indicators provide less dissipation and better resolution than the standard ones.
3.2.2. Smoothness indicators for nonuniform grids
Liu et al. (2018) defined the smoothness indicators (8) for uniform grids only. The generalization to the nonuniform case is presented in the following. For notational convenience
The firstorder threepoint numerical derivatives of y(x) at nodes x_{i − 1}, x_{i} and x_{i + 1} for nonuniform grids are approximated by (e.g., Singh & Bhadauria 2009)
Following the underlying idea of Liu et al. (2018), the new indicators are constructed as follows
If h_{i − 1} = h_{i}, the indicators (10) reduce to the uniform case indicators (8).
In monotonic smooth regions, the numerical derivatives (9) have the same sign, meaning that , and the indicators (10) reduce to
and, consequently,
In monotonic smooth regions, the nonlinear weights are exactly equal to the optimal weights, reducing numerical dissipation for the WENO interpolation.
Assume now that y(x) has a discontinuity in the interval [x_{i − 1}, x_{i}] and is smooth in the interval [x_{i}, x_{i + 1}]. In this case one has
Consequently, and
resulting in β_{1} ≫ β_{2}. Analogously, if y(x) is smooth in [x_{i − 1}, x_{i}] and has a discontinuity inside [x_{i}, x_{i + 1}], one obtains
resulting in β_{1} ≪ β_{2}. This guarantees that the smoothness indicators (10) effectively detect discontinuities in nonuniform grids and ensure the nonoscillatory property in the interpolation. The indicators (8) and (10) do not distinguish between discontinuities and smooth extrema. For this reason, they do not guarantee accuracy near such critical points. Some numerical evidence for the thirdorder WENO interpolation with smoothness indicators (10) is shown in Figs. 1–6.
Fig. 1. Exponential function (16), is approximated in the interval x ∈ [ − 1, 1] with thirdorder ENO, thirdorder WENO, cubic Lagrange, cubic Spline, monotonic cubic Hermite, and fourthorder WENO interpolations, with different homogeneously spaced grid points densities. Black dots represent the data values on the 36point grid. 
3.3. Fourthorder WENO interpolation
The WENO interpolation and the relative smoothness indicators outlined in the following are completely original. This novel fourthorder WENO interpolation uses a fourpoint stencil, is symmetric around the considered interval [x_{i}, x_{i + 1}], and provides a suitable compromise of the computation cost and the accuracy.
The unique cubic Lagrange polynomial p(x) interpolating the fourpoint large stencil
can be written as
that is, as a linear combination of the quadratic Lagrange interpolations
based, respectively, on the two threepoint substencils
The linear weights are given by
and satisfy γ_{2}(x)+γ_{3}(x)=1. The linear weights depend just on the grid geometry and not on the values of the function. The fourthorder WENO interpolation in the interval [x_{i}, x_{i + 1}] reads
where the nonlinear weights are defined as
with
and ϵ = 10^{−6}.
3.3.1. Smoothness indicators for uniform grids
The firstorder 4point numerical derivatives of y(x) at nodes x_{i − 1}, x_{i}, x_{i + 1} and x_{i + 2} for uniform grids read (e.g., Singh & Bhadauria 2009)
or, alternatively,
Using these, the new indicators are constructed as follows
with
In smooth regions where the second derivative has a constant sign, the differences of the numerical derivatives (13) have the same sign, that is,
Consequently, the indicators (12) reduce to
and one gets
Hence, in smooth regions where the second derivative has a constant sign, the nonlinear weights correspond to the linear weights, reducing the numerical dissipation of the fourthorder WENO interpolation. We note that the indicators (12) allow for over or undershoots around smooth extrema, yielding a higher accuracy near such critical points with respect to monotonic interpolations. However, this feature may have the adverse effect of producing negative interpolated values in strictly positive functions. Moreover, the indicators (12) do not guarantee accuracy near inflection points.
Order of accuracy.
For the discontinuous case, assume that y(x) has a discontinuity in the interval [x_{i − 1}, x_{i}] and is smooth in the interval [x_{i}, x_{i + 2}]. In this case one has
and from Eq. (13), one has
Consequently,
resulting in β_{2} ≫ β_{3}. Analogously, if the function y(x) is smooth in [x_{i − 1}, x_{i + 1}] and has a discontinuity inside [x_{i + 1}, x_{i + 2}], one obtains
resulting in β_{2} ≪ β_{3}. If the discontinuity point is located in the interval [x_{i}, x_{i + 1}], both substencils and contain the discontinuity. This seemingly difficult case is actually not problematic, because both polynomials q_{2}(x) and q_{3}(x) are essentially monotone inside [x_{i}, x_{i + 1}], because the over or undershoots appear in the cells adjacent to the discontinuity (Shu 1998, 2009).
3.3.2. Smoothness indicators for nonuniform grids
The firstorder fourpoint derivatives of y(x) at nodes x_{i − 1}, x_{i}, x_{i + 1} and x_{i + 2} for nonuniform grids are approximated by (e.g., Singh & Bhadauria 2009)^{5}
where H = h_{i − 1} + h_{i} + h_{i + 1}. The novel smoothness indicators are constructed as follows
For h_{i − 1} = h_{i} = h_{i + 1} = h, the finite difference formulae (14) and (15) reduce to uniform grids formulae (11) and (12), respectively.
In smooth regions where the second derivative has a constant sign, the differences of two consecutive numerical derivatives given by Eq. (14) have the same sign, that is,
In this case, one shows with some tedious algebra that the indicators (15) satisfy β_{2} = β_{3} and, consequently,
Hence, in smooth regions where the second derivative has a constant sign, the nonlinear weights correspond to the optimal weights, reducing numerical dissipation. As above, indicators (15) allow for over or undershoots around smooth extrema, yielding a higher accuracy near such critical points with respect to monotonic interpolations. As for the uniform case, this feature may have the adverse effect of producing negative interpolated values in strictly positive function.
The capability of the smoothness indicators (15) to detect discontinuities in the interpolation is already demonstrated for the uniform case in Sect. 3.3.1. The generalization of the proof to nonuniform grids is particularly cumbersome, but leads to the same conclusion. However, the indicators (15) do not guarantee accuracy near inflection points.
In conclusion, the smoothness indicators (12) and (15) effectively detect discontinuities and ensure the nonoscillatory property in the interpolation. Some numerical evidence for the fourthorder WENO interpolations with smoothness indicators (12) and (15) is displayed in Figs. 1–6. Moreover, an alternative fourthorder WENO interpolation is presented in Appendix A.
4. Numerical tests
For the sake of comparison, the thirdorder ENO, the thirdorder WENO (version Liu et al. 2018), the cubic Lagrange, the cubic Spline, the local monotonic piecewise cubic Hermite (Fritsch & Butland 1984; Ibgui et al. 2013), and the novel fourthorder WENO interpolations are tested. Four different functions are interpolated on both uniform and nonuniform grids with different numbers of grid points. The nonuniform grids are randomly generated. The experimental orders of convergence on uniform grids are summarized in Table 1. We note that the monotonic cubic Hermite interpolation is thirdorder accurate only. This is due to the derivatives provided by Fritsch & Butland (1984), which are secondorder accurate on uniform grids. In order to verify the accuracy for smooth cases, in Figs. 1 and 2 we analyze the interpolation of the exponential function
No significant difference between the various interpolations is visible for the homogeneous samplings, apart from the ENO technique that shows a slightly fragmented interpolation. Major differences appear for the inhomogeneous samplings, where fourthorder interpolations perform definitely better and the ENO technique shows a fragmented interpolation. Moreover, the monotonic cubic Hermite interpolation proves to be less accurate for inhomogeneous samplings. This is due to the derivatives provided by Fritsch & Butland (1984), which drop to firstorder accuracy on nonuniform grids.
In order to understand the behaviors around discontinuities, Figs. 3 and 4 show the interpolation of the scaled Heaviside function, that is,
The cubic Lagrange and the cubic Spline interpolations reveal Gibbslike oscillations for homogeneous and inhomogeneous samplings. For the homogeneous sampling, such oscillations do not decay in magnitude when the computational grid is refined. By contrast, monotonic cubic Hermite, ENO and WENO interpolations approximate the discontinuity more sharply and without oscillations. In this example, the WENO interpolation algorithm successfully capitalizes on its main idea, that is, placing much heavier weights on the candidate substencils in which the discontinuity is absent.
In order to better understand the different behaviors around smooth critical points (smooth extrema), Figs. 5 and 6 analyze the interpolation of the function
With 16 or more uniformly distributed grid points, thirdorder ENO, the cubic Lagrange, the cubic Spline interpolations, and WENO 4 accurately represent the minimum. All these methods require 36 points to accurately represent the extremum for inhomogeneous samplings. By contrast, the accuracy of WENO 3 and of the monotonic cubic Hermite interpolations locally decreases near the smooth extremum of the interpolated function. We note that the WENO 4 method allows for over or undershoots around the smooth extrema, yielding a higher accuracy with respect to monotonic interpolants near such critical points.
Figures 7 and 8 show the interpolation of the function
which contains both a highorder smooth structure and a discontinuity. Standard interpolations present oscillations for homogeneous and inhomogeneous samplings. The ENO interpolation is robust and resolves the complex discontinuous function quite well, albeit with small overshoots. WENO 3 and the monotonic cubic Hermite interpolation do not present oscillations, but struggle with accurately reproducing the smooth extrema in the function. By contrast, WENO 4 shows no oscillations, nor does it produce the spread of the extrema.
5. Application to radiative transfer
High accuracy is required in many radiative transfer problems (SocasNavarro et al. 2000; Trujillo Bueno 2003). In particular, the upcoming American Daniel K. Inouye Solar Telescope (DKIST; Elmore et al. 2014; Tritschler et al. 2016) and the planned European Solar Telescope (EST; Matthews et al. 2016; Matthews & Collados 2017) have a diffraction limit near to or surpassing the resolution of the best current 3D radiativemagnetohydrodynamic (RMHD) simulations of the solar surface (Peck et al. 2017), making high accuracy critical for comparisons with observations. Such 3D RMHD atmospheric models can be highly inhomogeneous and dynamic and present discontinuities. However, common radiative transfer calculations usually assume smooth variations in the radiation field and in the atmospheric physical parameters, and the incidence of discontinuities is often neglected (Steiner et al. 2016; Janett 2019). In fact, standard numerical schemes and highorder interpolations tend to misrepresent the nonsmooth regions of a problem, introducing spurious oscillations near discontinuities (e.g., Richards 1991; Zhang & Martin 1997; Shu 1998). Moreover, many radiative transfer codes require grid refinements, that is, the process of resolving the input data on a finer grid.
The fourthorder WENO interpolation technique presented in Sects. 3.3 is particularly suitable for scenarios that involve both complex smooth structures and discontinuities in the physical parameters. Moreover, it uses the same fourpoint stencil as the cubic Lagrange and monotonic cubic interpolations. Therefore, it is fairly straightforward to implement them on already existing codes.
5.1. 2D and 3D problems
The only difference between 1D and multidimensional radiative transfer lies in the mapping of the relevant quantities along the path of the photons (Carlsson 2008), typically through interpolation. At least for the Cartesian type grids, multidimensional interpolations can be obtained from 1D procedures. In 2D problems, the usual strategy is to discretize the ray by taking its intersections with the segments connecting the grid points, dealing then with 1D interpolations to recover offgrid points quantities (Auer 2003). In 3D problems, an efficient approach is to discretize the ray by taking its intersections with the individual planes defined by the grid points, dealing then with 2D interpolations.
Many radiative transfer codes have the option to use linear or bilinear interpolations from the values of the nearest grid points: for example, RH by Uitenbroek (2001), MULTI3D by Leenaarts & Carlsson (2009), and PORTA by Štěpán & Trujillo Bueno (2013). In multiple dimensions, the commonly used shortcharacteristic method suffers large numerical diffusion if the upwind intensity is successively linearly interpolated in the direction of propagation. Consequently, the solution lacks highorder accuracy (Fabiani Bendicho 2003; Ibgui et al. 2013; Peck et al. 2017). Such diffusion errors decrease with increasing order of accuracy of the interpolation and a natural choice would be to use biquadratic or bicubic interpolations (Kunasz & Auer 1988; Dullemond 2013)^{6}. However, if the behavior of the quantities to be interpolated is discontinuous or particularly intermittent, highorder interpolants introduce spurious extrema that could lead to nonphysical results (e.g., negative values of the source function or of the intensity). To avoid this, one can opt for monotonic interpolation schemes (e.g., Steffen 1990; Auer & Paletou 1994; Hayek et al. 2010; Ibgui et al. 2013), which however degenerate to a linear interpolation near smooth extrema, dropping to secondorder accuracy. This calls for highorder interpolations that are able to handle discontinuities and the WENO interpolation presented in Sect. 3.3 can be easily generalized to 2D and 3D Cartesian grids with a direct application dimension by dimension (Zhang et al. 2016).
5.2. Multigrid methods
Steiner (1991) first implemented a linear multigrid method for radiative transfer problems. Later on, Fabiani Bendicho et al. (1997) applied the nonlinear multigrid method to the multilevel NLTE radiative transfer problem. Recently, Štěpán & Trujillo Bueno (2013) generalized the method to the polarized case and Bjørgen & Leenaarts (2017) applied a nonlinear multigrid method to realistic 3D RMHD atmospheric models, which are highly inhomogeneous and dynamic.
The Jacobi and GaussSeidel iterative schemes for 3D NLTE radiative transfer problems quickly smooth out the highspatialfrequency error, whereas they slowly decrease the lowfrequency error. To accelerate the global smoothing of errors, multigrid methods transfer the problem to a coarser grid, where such lowfrequency errors become highfrequency errors and iterations on the coarse grid quickly decrease these errors. The coarser grid correction is mapped to the finer grid with an interpolation operator. In this step, cubic interpolations appear to give higher convergence rates (Štěpán & Trujillo Bueno 2013). However, when abrupt changes of atmospheric physical quantities are present (e.g., in 3D RMHD atmospheric models), such interpolations introduce spurious oscillations, which negatively affect the convergence rate of multigrid methods and even induce numerical instabilities, leading to unphysical negative values of strictly positive physical parameters. In such situations, it is safer to use monotonic interpolation, such as the monotonic cubic Hermite interpolation used in MULTI3D (Bjørgen & Leenaarts 2017). The interpolation operation is thus a crucial ingredient of any multigrid method and the fourthorder WENO interpolation might be an improvement over currently used interpolations.
5.3. Conversion to optical depth
The steadystate version of the radiative transfer equation consists in the linear firstorder inhomogeneous ODE given by (e.g., Mihalas 1978)
where the spatial coordinate s denotes the position along the ray under consideration, ν is the frequency, I_{ν} is the specific intensity, and χ_{ν} and ϵ_{ν} are the absorption and the emission coefficients, respectively. In nonlocal thermodynamic equilibrium (NLTE) conditions, the absorption and the emission coefficients depend in a complicated manner on the intensity I_{ν}, so that Eq. (20) is nonlinear and must be supplemented by additional statistical equilibrium equations and/or by suitable redistribution matrices for scattering processes.
The change of coordinates defined by
transplants Eq. (20) to the optical depth scale, namely,
where S_{ν} = ϵ_{ν}/χ_{ν} is the source function. From Eq. (21), one has
and numerical approximations of τ_{ν}(s) can be obtained by replacing the integral with a numerical quadrature. Janett et al. (2017a) explain that highorder formal solvers require a corresponding highorder quadrature of the integral above.
A npoint weighted quadrature rule is usually stated as
and is based on a suitable choice of the nodes {x_{i}} and the weights {ω_{i}}. In highorder quadratures, one must recover the values of f(x) at offgrid points, which typically requires highorder interpolations.
As a practical example, the fourthorder accurate Simpson’s formula
considers the node points and their corresponding function values. The midpoint function value must be recovered through a highorder interpolation, for example, the fourthorder WENO interpolation.
5.4. Highorder formal solvers
When applied to Eq. (20) or Eq. (22), certain highorder formal solvers, such as the RK4 method proposed by Landi Degl’Innocenti (1976) and the thirdorder pragmatic method by Janett et al. (2018), require the absorption χ_{ν} and emission ϵ_{ν} coefficients at offgrid points. In order to maintain highorder accuracy, such quantities must be obtained through highorder interpolations (e.g., Janett et al. 2017b, for the polarized problem), which are notoriously oscillatory in the presence of abrupt changes or discontinuities in the atmospheric physical quantities. In such formal solvers, spurious oscillations in the interpolation of χ_{ν} and ϵ_{ν} can negatively affect their order of accuracy. The fourthorder WENO interpolation might be befitting in such situations.
5.5. Redistribution matrix formalism
The redistribution matrix formalism is often used for calculating the emissivity in NLTE conditions, especially when partial frequency redistribution effects are important. The emissivity is evaluated by integrating the incident radiation field, multiplied by the redistribution matrix, over frequency and direction. However, performing the frequency integration is often nontrivial because of the sharp frequency dependence of the redistribution matrix. The quadrature of such integral requires describing the frequency dependence of the radiation field between the specific frequency points where its values are known. Interpolation techniques such as cardinal natural cubic splines have already been used (see Adams et al. 1971; Belluzzi & Trujillo Bueno 2014; Alsina Ballester et al. 2017). However, such a choice is problematic if the radiation field is not sufficiently smooth in the considered frequency interval. This calls for highorder interpolations capable of dealing with such sharp variations, such as the fourthorder WENO interpolation.
6. Conclusions
This paper discusses different ENO and WENO interpolation techniques. These interpolations are particularly suitable to problems containing both sharp discontinuities and complex smooth structures, where highorder Lagrange or spline interpolations are prone to over and undershoots.
Special attention is paid to the possible applications of a novel fourthorder WENO interpolation, which is outlined, analyzed, and tested for both the uniform and nonuniform cases. This method for computing nonoscillatory fourthorder interpolants is simple, symmetrical, and completely local. Unlike Bézier and monotonic Hermite interpolations, it does not degenerate to secondorder accuracy near smooth extrema of the interpolated function and it avoids the use of conditional “if…else” statements in the algorithm. Moreover, it uses the same fourpoint stencil as the cubic Lagrange interpolation. The implementation on already existing codes is, therefore, fairly straightforward. Numerical analysis and numerical experiments indicate that the smoothness indicators given by Eqs. (12) and (15) guarantee fourthorder accuracy in smooth regions, while effectively detecting discontinuities and enabling the nonoscillatory property. These novel local smoothness indicators yield less dissipation and better resolution than the canonical ones. Such indicators allow for minimal over or undershoots around smooth extrema, yielding a higher accuracy near such critical points. However, we note that this feature may have the adverse effect of producing negative interpolated values in strictly positive function. Moreover, this WENO interpolation is not differentiable at the grid points (i.e., it has kinks). Even though the fourthorder WENO interpolation may require more CPU time than the cubic Lagrange interpolation (depending on the specific algorithm components and computer implementation), it may still be computationally advantageous for many problems because of its high accuracy. Moreover, the method can be extended to the interpolation of 2D or 3D data.
This interpolation technique might be particularly suitable in the context of the numerical radiative transfer, which must be performed with ever increasing spatial and spectral resolution and consequently requires high accuracy in calculations. However, the fourthorder WENO interpolation is a general approximation procedure, which can also be used in many other applications beyond radiative transfer, such as computer vision and image processing. Alternative approximation techniques are available, such as the socalled HWENO approximations – WENO schemes based on Hermite polynomials (Qiu & Shu 2004; Zhu & Qiu 2009).
The formulae in Singh & Bhadauria (2009) contain some typos, which have been corrected in the present manuscript.
Acknowledgments
The financial support by the Swiss National Science Foundation (SNSF) through grant ID 200021_159206 is gratefully acknowledged. Special thanks are extended to A. Paganini and A. Battaglia for particularly enriching discussions. The authors are particularly grateful to the anonymous referee for providing valuable comments that helped improving the article.
References
 Adams, T. F., Hummer, D. G., & Rybicki, G. B. 1971, J. Quant. Spectr. Rad. Transf., 11, 1365 [NASA ADS] [CrossRef] [Google Scholar]
 Alsina Ballester, E., Belluzzi, L., & Trujillo Bueno, J. 2017, ApJ, 836, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Aràndiga, F., Baeza, A., & Yáñez, D. F. 2013, Adv. Compu. Math., 39, 289 [CrossRef] [Google Scholar]
 Auer, L. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 3 [NASA ADS] [Google Scholar]
 Auer, L. H., & Paletou, F. 1994, A&A, 285, 675 [NASA ADS] [Google Scholar]
 Belluzzi, L., & Trujillo Bueno, J. 2014, A&A, 564, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Bjørgen, J. P., & Leenaarts, J. 2017, A&A, 599, A118 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Borges, R., Carmona, M., Costa, B., & Don, W. S. 2008, J. Comput. Phys., 227, 3191 [Google Scholar]
 Carlsson, M. 2008, Phys. Scr. Vol. T, 133, 014012 [NASA ADS] [CrossRef] [Google Scholar]
 ČrnjarićŽic, N., Maćešić, S., & Crnković, B. 2007, Ann. Univ. Ferrara, 53, 199 [CrossRef] [Google Scholar]
 Dullemond, C. P. 2013, Radiative Transfer in Astrophysics. Theory, Numerical Methods and Applications, Tech. rep. (University of Heidelberg) [Google Scholar]
 Elmore, D. F., Rimmele, T., Casini, R., et al. 2014, in Groundbased and Airborne Instrumentation for Astronomy V, Proc. SPIE, 9147, 914707 [Google Scholar]
 Fabiani Bendicho, P. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 419 [NASA ADS] [Google Scholar]
 Fabiani Bendicho, P., Trujillo Bueno, J., & Auer, L. 1997, A&A, 324, 161 [Google Scholar]
 Fjordholm, U. S. 2016, in Handbook of Numerical Methods for Hyperbolic Problems: Basic and Fundamental Issues (Elsevier), Handbook of Numerical Analysis, 17, 123 [CrossRef] [Google Scholar]
 Fjordholm, U. S., Mishra, S., & Tadmor, E. 2011, ArXiv eprints [arXiv: 1112.1131] [Google Scholar]
 Fjordholm, U. S., Mishra, S., & Tadmor, E. 2012, SIAM J. Numer. Anal., 50, 544 [CrossRef] [Google Scholar]
 Fritsch, F. N., & Butland, J. 1984, SIAM J. Sci. Stat. Comput., 5, 300 [Google Scholar]
 Fritsch, F. N., & Carlson, R. E. 1980, SIAM J. Numer. Anal., 17, 238 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [Google Scholar]
 Gande, N. R., Rathod, Y., & Rathan, S. 2017, Int. J. Numer. Methods Fluids, 85, 90 [NASA ADS] [CrossRef] [Google Scholar]
 Harten, A., Osher, S., Engquist, B., & Chakravarthy, S. R. 1986, Appl. Numer. Math., 2, 347 [CrossRef] [MathSciNet] [Google Scholar]
 Harten, A., Engquist, B., Osher, S., & Chakravarthy, S. R. 1987, J. Comput. Phys., 71, 231 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Hayek, W., Asplund, M., Carlsson, M., et al. 2010, A&A, 517, A49 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Henrick, A. K., Aslam, T. D., & Powers, J. M. 2005, J. Comput. Phys., 207, 542 [NASA ADS] [CrossRef] [Google Scholar]
 Ibgui, L., Hubeny, I., Lanz, T., & Stehlé, C. 2013, A&A, 549, A126 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Janett, G. 2019, A&A, 622, A162 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Janett, G., Carlin, E. S., Steiner, O., & Belluzzi, L. 2017a, ApJ, 840, 107 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Steiner, O., & Belluzzi, L. 2017b, ApJ, 845, 104 [NASA ADS] [CrossRef] [Google Scholar]
 Janett, G., Steiner, O., & Belluzzi, L. 2018, ApJ, 865, 16 [NASA ADS] [CrossRef] [Google Scholar]
 Jiang, G.S., & Shu, C.W. 1996, J. Comput. Phys., 126, 202 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Kunasz, P., & Auer, L. H. 1988, J. Quant. Spectr. Rad. Transf., 39, 67 [Google Scholar]
 Landi Degl’Innocenti, E. 1976, A&A, 25, 379 [Google Scholar]
 Leenaarts, J., & Carlsson, M. 2009, in The Second Hinode Science Meeting: Beyond DiscoveryToward Understanding, eds. B. Lites, M. Cheung, T. Magara, J. Mariska, & K. Reeves, ASP Conf. Ser., 415, 87 [NASA ADS] [Google Scholar]
 Liu, S., Shen, Y., Chen, B., & Zeng, F. 2018, Int. J. Numer. Methods Fluids, 87, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, X.D., Osher, S., & Chan, T. 1994, J. Comput. Phys., 115, 200 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Matthews, S. A., & Collados, M. 2017, SOLARNET IV: The Physics of the Sun from the Interior to the Outer Atmosphere, 78 [Google Scholar]
 Matthews, S. A., Collados, M., Mathioudakis, M., & Erdelyi, R. 2016, in Groundbased and Airborne Instrumentation for Astronomy VI, Proc. SPIE, 9908, 990809 [Google Scholar]
 Mihalas, D. 1978, Stellar Atmospheres, 2nd edn. (San Francisco: W. H. Freeman and Company) [Google Scholar]
 Peck, C. L., Criscuoli, S., & Rast, M. P. 2017, ApJ, 850, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Qiu, J., & Shu, C.W. 2004, J. Comput. Phys., 193, 115 [NASA ADS] [CrossRef] [Google Scholar]
 Richards, F. 1991, J. Approximation Theory, 66, 334 [CrossRef] [Google Scholar]
 Shu, C. W. 1998, Essentially nonoscillatory and weighted essentially nonoscillatory schemes for hyperbolic conservation laws (Springer), 325 [Google Scholar]
 Shu, C.W. 2009, SIAM Rev., 51, 82 [NASA ADS] [CrossRef] [Google Scholar]
 Shu, C.W., & Osher, S. 1988, J. Comput. Phys., 77, 439 [Google Scholar]
 Shu, C.W., & Osher, S. 1989, J. Comput. Phys., 83, 32 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Singh, A. K., & Bhadauria, B. S. 2009, Int. J. Math. Anal., 3, 815 [Google Scholar]
 SocasNavarro, H., Bueno, J. T., & Cobo, B. R. 2000, ApJ, 530, 977 [NASA ADS] [CrossRef] [Google Scholar]
 Steffen, M. 1990, A&A, 239, 443 [NASA ADS] [Google Scholar]
 Steiner, O. 1991, A&A, 242, 290 [NASA ADS] [Google Scholar]
 Steiner, O., Züger, F., & Belluzzi, L. 2016, A&A, 586, A42 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Štěpán, J., & Trujillo Bueno, J. 2013, A&A, 557, A143 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Tritschler, A., Rimmele, T. R., Berukoff, S., et al. 2016, Astron. Nachr., 337, 1064 [NASA ADS] [CrossRef] [Google Scholar]
 Trujillo Bueno, J. 2003, in Stellar Atmosphere Modeling, eds. I. Hubeny, D. Mihalas, & K. Werner, ASP Conf. Ser., 288, 551 [NASA ADS] [Google Scholar]
 Uitenbroek, H. 2001, ApJ, 557, 389 [NASA ADS] [CrossRef] [Google Scholar]
 van Leer, B. 1979, J. Comput. Phys., 32, 101 [Google Scholar]
 Wu, X., Liang, J., & Zhao, Y. 2015, Int. J. Numer. Methods Fluids, 81 [Google Scholar]
 Zhang, Y. T., & Shu, C. W. 2016, in Handbook of Numerical Methods for Hyperbolic Problems, ed. R. Abgrall, & C. W. Shu (Elsevier), Handbook of Numerical Analysis, 17, 103 [Google Scholar]
 Zhang, Z., & Martin, C. F. 1997, J. Comput. Applied Math., 87, 359 [CrossRef] [Google Scholar]
 Zhu, J., & Qiu, J. 2009, J. Sci. Comput., 39, 293 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: An alternative version of the fourthorder WENO interpolation
In this appendix, an alternative to the fourthorder WENO interpolation presented in Sect. 3.3 is outlined. The unique cubic Lagrange polynomial p(x) interpolating the fourpoint large stencil
can be written as
in other words, as a linear combination of the linear Lagrange interpolations
based, respectively, on the three twopoint substencils
The linear weights are given by
and satisfy γ_{2}(x)+γ_{3}(x)+γ_{4}(x)=1. The linear weights depend just on the grid geometry and not on the function values. The fourthorder WENO interpolation in the interval [x_{i}, x_{i + 1}] reads
where the nonlinear weights are defined as
with
and ϵ = 10^{−6}.
A.1. Smoothness indicators for uniform grids
In case of uniform grids, a possible choice for the smoothness indicators is
where the differences of the numerical derivatives are given by Eq. (13). Unfortunately, this version makes the fourthorder WENO interpolation too dissipative and generates apparent oscillations. The search of smoothness indicators that yield less dissipation and better resolution is not over.
A.2. Smoothness indicators for nonuniform grids
In case of nonuniform grids, a possible choice for the smoothness indicators is
where the numerical derivatives are given by Eq. (14). As for the uniform case, this version makes the fourthorder WENO interpolation too dissipative and generates apparent oscillations. The search of smoothness indicators that yield less dissipation and better resolution is not over.
All Tables
All Figures
Fig. 1. Exponential function (16), is approximated in the interval x ∈ [ − 1, 1] with thirdorder ENO, thirdorder WENO, cubic Lagrange, cubic Spline, monotonic cubic Hermite, and fourthorder WENO interpolations, with different homogeneously spaced grid points densities. Black dots represent the data values on the 36point grid. 

In the text 
Fig. 2. Same as Fig. 1, but for inhomogeneous samplings. 

In the text 
Fig. 3. Same as Fig. 1, but for the modified Heaviside function given by Eq. (17). 

In the text 
Fig. 4. Same as Fig. 3, but for inhomogeneous samplings. 

In the text 
Fig. 5. Same as Fig. 1, but for the modified Gaussian function given by Eq. (18). 

In the text 
Fig. 6. Same as Fig. 5, but for inhomogeneous samplings. 

In the text 
Fig. 7. Same as Fig. 1, but for the discontinuous sine function given by Eq. (19). 

In the text 
Fig. 8. Same as Fig. 7, but for inhomogeneous samplings. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.