© 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 non-oscillatory fashion and at the same time to perform with uniform high-order accuracy in smooth regions.

It is a well known fact that standard high-order 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, monotonicity-preserving strategies such as Bézier and monotonic high-order 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 second-order accuracy (Shu 1998). By contrast, essentially non-oscillatory (ENO) and weighted ENO (WENO) techniques accomplish the feat of approximating a function with high-order accuracy in smooth parts, while avoiding “Gibbs-like” oscillations near the discontinuities (Fjordholm 2016).

Section 2 presents the ENO approximation technique and it exposes the third-order ENO interpolation. Section 3 presents the WENO approximation technique, with a special focus on a novel fourth-order WENO interpolation. Section 4 exposes numerical tests for different interpolations. Section 5 describes different radiative transfer problems that require high-order 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 (second-order accurate) to a piecewise constant approximation (first-order 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 first-order derivatives (Fritsch & Carlson 1980), whereas Bézier interpolants make use of the so-called 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 monotonicity-preserving 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², that is, twice continuously differentiable.

By contrast, ENO techniques are high-order approximations that handle discontinuities and maintain high-order 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 non-oscillatory property. There is broad and well-cited 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 j-th 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¹

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 p-th order ENO interpolation E(x) in the interval [x₁, x_N] reads

(1)

with

(2)

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 p-point stencil containing x_i is given by

(3)

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 one-point stencil

and adaptively adds points with the procedure described above, ending up with a p-point 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

(4)

where the Lagrange basis polynomials ℓ_k are defined as²

Hence, once the optimal candidate substencil is determined, the p-order 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 non-oscillatory 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 third-order accurate ENO interpolation is presented in the following.

2.3. Third-order ENO interpolation

In order to build a third-order accurate p_i(x), one needs a three-point stencil. Thus, starting from the one-point 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 two-point stencil is taken as

otherwise,

Suppose that the two-point 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 three-point stencil is taken as

otherwise,

Then, one proceeds with the Lagrange interpolation given by Eq. (4) on the selected substencil, yielding the third-order accurate quadratic polynomial p_i(x).

We note that the choice of the two-point stencil at the previous step may lead to the three-point stencil

In order to build p_i(x), the third-order 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 high-order accuracy in smooth regions while keeping the essentially non-oscillatory 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 arbitrary-order accurate WENO approximations for finite difference schemes and most applications use their fifth-order 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 p-th order WENO interpolation W(x) in the interval [x₁, x_N] reads

(5)

where I_i(x) is defined in Eq. (2). Consider the p-point stencil given by Eq. (3) to be the large stencil. The unique p-th-order accurate Lagrange polynomial p(x) interpolating can be recovered as a linear combination of ℓ-th-order accurate Lagrange polynomials q_m(x) defined on a set of ℓ-point substencils with ℓ ∈ {1, …, p} and m ∈ { − p + q + ℓ, …, q}, i.e.,³

(6)

where γ_m(x) are the so-called 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⁻⁶ (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)

(7)

is well established. Unfortunately, indicators (7) consider uniform grids alone and are particularly dissipative (i.e., they smear sharp gradients) for low-order WENO interpolations. In the following, a third-order and a novel fourth-order accurate WENO interpolations are presented.

3.2. Third-order WENO interpolation

Compared with higher-order versions, the third-order WENO interpolation is more robust for the treatment of discontinuous problems, it uses fewer grid points⁴, and it provides a suitable compromise between computational cost and accuracy (Liu et al. 2018).

The unique quadratic Lagrange polynomial p(x) interpolating the three-point large stencil

can be written as

in other words, as a linear combination of the linear Lagrange interpolations

based, respectively, on the two two-point substencils

The linear weights are given by

and satisfy γ₁(x)+γ₂(x)=1. The linear weights depend just on the grid geometry and not on the function values. The third-order WENO interpolation in the interval [x_i, x_i + 1] reads

where the nonlinear weights are defined as

with

and ϵ = 10⁻⁶.

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 third-order WENO interpolation too dissipative.

Alternatively, the mapping function of the WENO-M scheme (Henrick et al. 2005) and the global stencil indicator of the WENO-Z scheme (Borges et al. 2008) effectively improve the accuracy of high-order WENO schemes, but they are not satisfactory for the third-order one. More recently, Wu et al. (2015) proposed the global smoothness indicators for the less dissipative third-order WENO-N3 and WENO-NP3 schemes. Similarly, Gande et al. (2017) proposed the global smoothness indicator of the WENO-F3 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,

(8)

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 first-order three-point 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)

(9)

Following the underlying idea of Liu et al. (2018), the new indicators are constructed as follows

(10)

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 β₁ ≫ β₂. 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 β₁ ≪ β₂. This guarantees that the smoothness indicators (10) effectively detect discontinuities in nonuniform grids and ensure the non-oscillatory 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 third-order 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 third-order ENO, third-order WENO, cubic Lagrange, cubic Spline, monotonic cubic Hermite, and fourth-order WENO interpolations, with different homogeneously spaced grid points densities. Black dots represent the data values on the 36-point grid.

Fig. 2. Same as Fig. 1, but for inhomogeneous samplings.

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

Fig. 4. Same as Fig. 3, but for inhomogeneous samplings.

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

Fig. 6. Same as Fig. 5, but for inhomogeneous samplings.

3.3. Fourth-order WENO interpolation

The WENO interpolation and the relative smoothness indicators outlined in the following are completely original. This novel fourth-order WENO interpolation uses a four-point 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 four-point large stencil

can be written as

that is, as a linear combination of the quadratic Lagrange interpolations

based, respectively, on the two three-point substencils

The linear weights are given by

and satisfy γ₂(x)+γ₃(x)=1. The linear weights depend just on the grid geometry and not on the values of the function. The fourth-order WENO interpolation in the interval [x_i, x_i + 1] reads

where the nonlinear weights are defined as

with

and ϵ = 10⁻⁶.

3.3.1. Smoothness indicators for uniform grids

The first-order 4-point 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)

(11)

or, alternatively,

Using these, the new indicators are constructed as follows

(12)

with

(13)

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 fourth-order 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.

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 β₂ ≫ β₃. 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 β₂ ≪ β₃. 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₂(x) and q₃(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 first-order four-point 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)⁵

(14)

where H = h_i − 1 + h_i + h_i + 1. The novel smoothness indicators are constructed as follows

(15)

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 β₂ = β₃ 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 non-oscillatory property in the interpolation. Some numerical evidence for the fourth-order WENO interpolations with smoothness indicators (12) and (15) is displayed in Figs. 1–6. Moreover, an alternative fourth-order WENO interpolation is presented in Appendix A.

4. Numerical tests

For the sake of comparison, the third-order ENO, the third-order 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 fourth-order 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 third-order accurate only. This is due to the derivatives provided by Fritsch & Butland (1984), which are second-order 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

(16)

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 fourth-order 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 first-order accuracy on non-uniform grids.

In order to understand the behaviors around discontinuities, Figs. 3 and 4 show the interpolation of the scaled Heaviside function, that is,

(17)

The cubic Lagrange and the cubic Spline interpolations reveal Gibbs-like 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

(18)

With 16 or more uniformly distributed grid points, third-order 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

(19)

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

Fig. 8. Same as Fig. 7, but for inhomogeneous samplings.

which contains both a high-order 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 (Socas-Navarro 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 radiative-magnetohydrodynamic (R-MHD) simulations of the solar surface (Peck et al. 2017), making high accuracy critical for comparisons with observations. Such 3D R-MHD 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 high-order 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 fourth-order 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 four-point 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 off-grid 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 bi-linear 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 short-characteristic method suffers large numerical diffusion if the upwind intensity is successively linearly interpolated in the direction of propagation. Consequently, the solution lacks high-order 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 bi-quadratic or bi-cubic interpolations (Kunasz & Auer 1988; Dullemond 2013)⁶. However, if the behavior of the quantities to be interpolated is discontinuous or particularly intermittent, high-order 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 second-order accuracy. This calls for high-order 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 multi-level 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 R-MHD atmospheric models, which are highly inhomogeneous and dynamic.

The Jacobi and Gauss-Seidel iterative schemes for 3D NLTE radiative transfer problems quickly smooth out the high-spatial-frequency error, whereas they slowly decrease the low-frequency error. To accelerate the global smoothing of errors, multigrid methods transfer the problem to a coarser grid, where such low-frequency errors become high-frequency 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 R-MHD 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 fourth-order WENO interpolation might be an improvement over currently used interpolations.

5.3. Conversion to optical depth

The steady-state version of the radiative transfer equation consists in the linear first-order inhomogeneous ODE given by (e.g., Mihalas 1978)

(20)

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 non-local 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

(21)

transplants Eq. (20) to the optical depth scale, namely,

(22)

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 high-order formal solvers require a corresponding high-order quadrature of the integral above.

A n-point weighted quadrature rule is usually stated as

and is based on a suitable choice of the nodes {x_i} and the weights {ω_i}. In high-order quadratures, one must recover the values of f(x) at off-grid points, which typically requires high-order interpolations.

As a practical example, the fourth-order accurate Simpson’s formula

considers the node points and their corresponding function values. The midpoint function value must be recovered through a high-order interpolation, for example, the fourth-order WENO interpolation.

5.4. High-order formal solvers

When applied to Eq. (20) or Eq. (22), certain high-order formal solvers, such as the RK4 method proposed by Landi Degl’Innocenti (1976) and the third-order pragmatic method by Janett et al. (2018), require the absorption χ_ν and emission ϵ_ν coefficients at off-grid points. In order to maintain high-order accuracy, such quantities must be obtained through high-order 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 fourth-order 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 high-order interpolations capable of dealing with such sharp variations, such as the fourth-order 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 high-order Lagrange or spline interpolations are prone to over- and undershoots.

Special attention is paid to the possible applications of a novel fourth-order WENO interpolation, which is outlined, analyzed, and tested for both the uniform and nonuniform cases. This method for computing non-oscillatory fourth-order interpolants is simple, symmetrical, and completely local. Unlike Bézier and monotonic Hermite interpolations, it does not degenerate to second-order 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 four-point 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 fourth-order accuracy in smooth regions, while effectively detecting discontinuities and enabling the non-oscillatory 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 fourth-order 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 fourth-order 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 so-called HWENO approximations – WENO schemes based on Hermite polynomials (Qiu & Shu 2004; Zhu & Qiu 2009).

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

Appendix A: An alternative version of the fourth-order WENO interpolation

In this appendix, an alternative to the fourth-order WENO interpolation presented in Sect. 3.3 is outlined. The unique cubic Lagrange polynomial p(x) interpolating the four-point large stencil

can be written as

in other words, as a linear combination of the linear Lagrange interpolations

based, respectively, on the three two-point substencils

The linear weights are given by

and satisfy γ₂(x)+γ₃(x)+γ₄(x)=1. The linear weights depend just on the grid geometry and not on the function values. The fourth-order WENO interpolation in the interval [x_i, x_i + 1] reads

where the nonlinear weights are defined as

with

and ϵ = 10⁻⁶.

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 fourth-order 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 fourth-order 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 third-order ENO, third-order WENO, cubic Lagrange, cubic Spline, monotonic cubic Hermite, and fourth-order WENO interpolations, with different homogeneously spaced grid points densities. Black dots represent the data values on the 36-point 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