© ESO 2016

1. Introduction

Many stars, such as asymptotic giant branch (AGB) stars, have extended atmospheres, which has important physical and observational implications. In particular, the radiation field at large distances from the inner stellar disk becomes very dilute and is confined primarily to a narrow solid angle around the radial direction (Mihalas 1978). Additionally, planets can also feature extended atmospheres, in which case the use of the usual assumption of a plane-parallel atmosphere is no longer valid, for example, the extended and spherically stratified atmosphere of Pluto (Gladstone et al. 2016). Compared to the radiative transfer equation in the plane-parallel approximation, the transfer equation in spherical symmetry is more complicated and more difficult to handle. This is especially caused by the appearance of a derivative with respect to the polar angle, which is not present in the plane-parallel geometry. Other cases in which the plane-parallel radiative transfer equation becomes much more complicated include systems with a moving medium where photons are subject to the Doppler effect and aberration (Mihalas 1978).

Depending on the transport coefficients, the general radiative transport equation can have a very different mathematical character. As described in Kanschat et al. (2009), it is a hyperbolic wave equation in regions without matter, an elliptic diffusion equation in case of an optically thick medium, or a parabolic equation when the light is strongly peaked in the forward direction. It is, therefore, very difficult to find numerical methods which are capable of dealing with these different (or even mixed) types of behaviour.

The problem of a grey spherical atmosphere in radiative equilibrium was first studied by Kosirev (1934) and Chandrasekhar (1934). An extension of the iterative moment method for plane-parallel atmospheres using variable Eddington factors (Auer & Mihalas 1970) to the problem of a spherically symmetric atmosphere restricted to isotropic scattering was introduced by Hummer & Rybicki (1971). Additionally, solutions based on a generalised Eddington approximation were developed by Lucy (1971) and Unno & Kondo (1976), respectively.

Alternatively, the long (Cannon 1970) or short characteristic methods (Kunasz & Auer 1988) can be used. Combined with, for example, an operator perturbation method to account for the unknown source function – such as the accelerated lamdba iteration (Cannon 1973) – they can be used to solve the radiative transfer equation in multiple dimensions (Olson & Kunasz 1987; Kunasz & Olson 1988). In combination with the work of Hummer & Rybicki (1971), operator splitting methods have also been directly utilised for spherically symmetric atmospheres by Kubat (1994).

Currently, the computationally demanding Monte Carlo techniques can very efficiently solve the radiative transfer problem for a broad variety of applications ranging from stellar atmospheres (e.g. Abbott & Lucy 1985; Lucy & Abbott 1993; Lucy 1999), circumstellar disks (e.g. Pascucci et al. 2004, and the references therein), and supernovae (e.g. Lucy 2005) to dusty objects (e.g. Steinacker et al. 2013, and the references therein). Their advantages are, for example, that they are intrinsically three-dimensional, can handle easily complex geometries and density distributions, and have a low algorithmic complexity. However, in some special cases with very optically thick regions the Monte Carlo approach becomes very slow because the number of photon interactions (i.e. scattering) increases exponentially. In these situations grid-based solution techniques (e.g. a finite difference, finite element, or finite volume method) for the differential equation of radiative transfer are more favourable even if they are algorithmically quite complex. Moreover, grid-based techniques allow for higher order approximations as well as error control, and they can be formulated to be intrinsically flux conserving (Hesthaven & Warburton 2008; Henning 2001).

In this study, we present a numerical method for a direct solution of the radiative transfer equation in spherical symmetry by means of a finite element method. Finite element methods have already been used to solve the radiative transfer equation, most notably in the three-dimensional case by Kanschat (1996) or Richling et al. (2001). In these cases, a continuous Galerkin finite element was used to discretise the transport operator. Such an approach, however, causes problems in cases where the transport equation has a dominant hyperbolic character, which may lead to discontinuities within the solution. It is possible to stabilise the numerical scheme by introducing a streamline diffusion method, for example, which essentially transforms the hyperbolic equation into a parabolic one by adding a small diffusion in the direction of photon transport. While this stabilises the numerical discretisation, it also adds a free parameter to the method, which has to be chosen carefully in such a way that it stabilises the numerical scheme, but does not result in a decreased accuracy of the solution (see e.g. Eriksson et al. 1996; Kanschat 1996). An alternative way to deal with such problems is to soften the requirements for the solution of the finite element method, in particular the requirement of continuity of the solution across element boundaries. One such numerical approach in the context of finite element methods is the discontinuous Galerkin method.

The discontinuous Galerkin finite element method (DG-FEM) was first introduced by Reed & Hill (1973) to study neutron transport problems. Since then, the method has been successfully applied to a broad variety of parabolic, hyperbolic, and elliptic problems (e.g. Hesthaven & Warburton 2008; Cockburn 2003). For three-dimensional internal heat transfer problems, the discretisation with a discontinuous Galerkin method has been studied by Cui & Li (2004). For astrophysical applications, discontinuous Galerkin methods have, for example, been used by Dykema et al. (1996) and Castor et al. (1992).

In this work we use a DG-FEM to directly solve the radiative transfer equation with arbitrarily complex scattering phase functions in spherical symmetry as a two-dimensional problem. In contrast to other solution methods, the scattering integral can thereby be evaluated exactly, independent of the angular mesh resolution. The development of the numerical scheme is outlined for the formulation of the DG approach in Sects. 2 and 3. Its applicability to several test problems is presented in Sect. 4. Remarks on the numerical efficiency are given in Sect. 5, followed by an outlook (Sect. 6) and a summary (Sect. 7).

2. Radiative transfer equation

We consider the radiative transfer equation in spherical symmetry, cf. Mihalas (1978), i.e. $\begin{array}{ccc}\mathit{\mu }\hspace{0.17em}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,\mu }\mathrm{\right)}}{\mathit{\partial r}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,\mu }\mathrm{\right)}}{\mathit{\partial \mu }}\mathrm{=}\mathrm{-}{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,\mu }\mathrm{\right)}\mathrm{-}{\overline{\mathit{s}}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,\mu }\mathrm{\right)}& & \\ \mathrm{+}{\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{ind}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,\mu }\mathrm{\right)}\mathrm{+}{\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{sp}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& & \\ \mathrm{+}\frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{s}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{;}\mathit{\phi ,}{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}{\mathit{\phi }}^{\mathrm{\prime }}& & \end{array}$(1)with the specific intensity I_ν, radius coordinate r, the cosine of the polar angle μ, the azimuth angle φ, absorption coefficient κ_ν, total (${\overline{\mathit{s}}}_{\mathit{\nu }}$) and differential (s_ν(r,μ,μ′;φ,φ′)) scattering coefficients, and induced and spontaneous emission coefficients ${\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{ind}}\mathit{,}{\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{sp}}$.

The differential scattering coefficient can be decomposed in an averaged angle-independent part, the total scattering coefficient ${\overline{\mathit{s}}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r}\mathrm{\right)}$, and the phase function p(r,μ,μ′;φ,φ′). This yields $\begin{array}{ccc}{\mathit{s}}_{\mathit{\nu }}& & \mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{;}\mathit{\phi ,}{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\\ & & \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}{\mathit{s}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{;}\mathit{\phi ,}{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\hspace{0.17em}\mathrm{d}{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{·}\mathit{p}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{;}\mathit{\phi ,}{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{\right)}\\ & & \mathrm{=}{\overline{\mathit{s}}}_{\mathit{\nu }}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{·}\mathit{p}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{;}\mathit{\phi ,}{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{.}\end{array}$(2)The phase function describes the probability of photon scattering of incident directions μ′ and φ′ into the directions μ and φ. Since in spherical symmetry the radiation field is independent of the azimuth angle φ, the phase function p(r,μ,μ′;φ,φ′) can formally be integrated with respect to the azimuth angle to obtain the azimuthally averaged phase function ${\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathit{,\mu }\mathrm{\right)}\mathrm{:}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\pi }}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathit{p}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathit{,\mu }\mathrm{;}\mathit{\phi ,}{\mathit{\phi }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\phi }}^{\mathrm{\prime }}\mathit{,}$(3)which is then used in the transfer equation. To guarantee the conservation of energy, the phase function is normalised to unity, i.e. $\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{=}\mathrm{1.}$(4)The extinction coefficient ${\mathit{\chi }}_{\mathit{\nu }}\mathrm{=}{\mathit{\kappa }}_{\mathit{\nu }}\mathrm{+}{\overline{\mathit{s}}}_{\mathit{\nu }}$ can be combined with the induced emission into $\mathit{\chi ̂}\mathit{\nu }\mathrm{=}{\mathit{\chi }}_{\mathit{\nu }}\mathrm{-}{\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{ind}}\mathit{.}$(5)In this work we only consider cases with ${\mathit{\chi }}_{\mathit{\nu }}\mathit{>}{\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{ind}}$ and therefore , i.e. we do not consider systems which are dominated by induced emission (such as maser). Equation (1)is therefore simplified to $\mathit{\mu }\hspace{0.17em}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\nu }}}{\mathit{\partial r}}\mathrm{+}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{\nu }}}{\mathit{\partial \mu }}\mathrm{=}\mathrm{-}\mathit{\chi ̂}\mathit{\nu }\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\mathrm{+}{\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{sp}}\mathrm{+}\frac{{\overline{\mathit{s}}}_{\mathit{\nu }}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\mathit{,}$(6)where all direct dependences on r and μ have been omitted to simplify the presentation. Instead of the usual form of the radiative transfer equation in spherical symmetry given by Eq. (6), we use the mathematically equivalent form $\begin{array}{ccc}\frac{\mathit{\partial }}{\mathit{\partial r}}\left(\mathit{\mu }\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\right)\mathrm{+}\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}{\mathit{I}}_{\mathit{\nu }}\mathrm{=}& & \\ \mathrm{-}\mathit{\chi ̂}\mathit{\nu }\hspace{0.17em}{\mathit{I}}_{\mathit{\nu }}\mathrm{+}{\mathit{\eta }}_{\mathit{\nu }}^{\mathrm{sp}}\mathrm{+}\frac{{\overline{\mathit{s}}}_{\mathit{\nu }}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{I}}_{\mathit{\nu }}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\mathit{.}& & \end{array}$(7)To simplify the notation we drop the explicit frequency dependence in the transfer equation and its corresponding coefficients in the following. We do not introduce the usual (radial) optical depth in the transfer equation. Doing so – by moving the extinction coefficient to the left-hand side – would introduce additional terms of the form , which can be numerically destabilising since the extinction coefficient can vary over several orders of magnitude.

Fig. 1 Computational domain Ω. The boundaries of the domain between the inner radius rin and outer radius rout, as well as μ between + 1 and −1, are denoted by the respective Γ.

The Eq. (7)is solved on the two-dimensional domain Ω ⊂ R² with Ω ∋ x = (r,μ) ∈ (r_in,r_out) × (−1,1) as displayed in Fig. 1. The inner (outer) radius of Ω is given by r_in (r_out).

At the inner boundary Γ_in of the domain boundary conditions need to be imposed for positive values of μ, i.e. $\mathit{I}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{in}} \mathrm{for}{x}\mathrm{\in }{\mathrm{\Gamma }}_{\mathrm{in}}^{\mathrm{+}}$(8)given by the radiation field, for example, of the stellar interior. Likewise, at the outer boundary Γ_out an incident external radiation field has to be prescribed for negative values of μ, i.e. $\mathit{I}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{\mathit{I}}_{\mathrm{out}} \mathrm{for}{x}\mathrm{\in }{\mathrm{\Gamma }}_{\mathrm{out}}^{\mathrm{-}}\mathit{.}$(9)We note that since the differential with respect to μ in Eq. (7)vanishes on the boundary ∂Γ, no explicit boundary conditions need to be imposed on Γ_μ^±. Therefore, the boundary value problem considered in this work is given by $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{\mathit{\partial }}{\mathit{\partial r}}\left(\mathit{\mu }\hspace{0.17em}\mathit{I}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}\mathit{I}\right)\mathrm{=}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{-}\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\hspace{0.17em}\mathit{I}\mathrm{-}\mathit{\chi ̂}\hspace{0.17em}\mathit{I}\mathrm{+}{\mathit{\eta }}^{\mathrm{sp}}\mathrm{+}\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathit{I}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\prime \begin{array}{c}\end{array}\mathrm{in}\mathrm{\Omega }\mathit{,}\\ \begin{array}{c}\end{array}\frac{\mathit{\partial I}}{\mathit{\partial r}}\mathrm{=}\mathrm{-}\frac{\mathrm{2}}{\mathit{r}}\hspace{0.17em}\mathit{I}\mathrm{\mp }\mathit{\chi ̂}\hspace{0.17em}\mathit{I}\mathrm{\pm }{\mathit{\eta }}^{\mathrm{sp}}\mathrm{\pm }\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathit{I}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{(}\mathrm{\pm }\mathrm{1}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\prime \begin{array}{c}\end{array}\mathrm{on}{\mathrm{\Gamma }}_{{\mathit{\mu }}^{\mathrm{\pm }}}\mathit{,}\\ \begin{array}{c}\end{array}\mathit{I}\mathrm{=}{\mathit{I}}_{\mathrm{in}}\begin{array}{c}\end{array}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{in}}^{\mathrm{+}}\mathit{,}\\ \begin{array}{c}\end{array}\mathit{I}\mathrm{=}{\mathit{I}}_{\mathrm{out}}\begin{array}{c}\end{array}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{out}}^{\mathrm{-}}\mathit{.}\end{array}$(10)

3. Formulation of the discontinuous Galerkin method

For the formulation of the discontinuous Galerkin method we follow the notation of Hesthaven & Warburton (2008). The domain Ω is described by K non-overlapping elements D^k, i.e. $\mathrm{int}\mathrm{\left(}{\mathrm{D}}^{\mathit{k}}\mathrm{\right)}\mathrm{\cap }\mathrm{int}\mathrm{\left(}{\mathrm{D}}^{\mathit{n}}\mathrm{\right)}\mathrm{=}\mathrm{\varnothing } \mathrm{for}\mathit{k}\ne \mathit{n,}$(11)such that the continuous domain Ω is approximated by the finite domain Ω_h with $\mathrm{\Omega }\backsimeq {\mathrm{\Omega }}_{\mathit{h}}\mathrm{=}\underset{\mathit{k}\mathrm{=}\mathrm{1}}{\overset{\mathit{K}}{\mathrm{􏽛}}}{\mathrm{D}}^{\mathit{k}}\mathrm{=}\overline{\mathrm{\Omega }}$(12)with a local element size h^k: = diam(D^k).

Definition 1. For any open set Ω in Rⁿ let P_q(Ω) and Q_q(Ω) denote the set of algebraic polynomials of total and partial degree ≤q, respectively. For n = 2 Table 1 shows the monomial basis functions of P_q and Q_q for q = 0,1,2,3 and Fig. 2 displays their Lagrangian elements.

Fig. 2 Lagrangian reference elements for dimension n = 2, and polynomial degrees q = 0,1,2,3. The degrees of freedom (see Table 1) are marked with •. Top: Pq elements, bottom: Qq elements.

In this work we use the nodal representation of the solution ${\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{:}\mathrm{=}{\mathit{I}}_{\mathit{h}}{\mathrm{|}}_{{\mathrm{D}}^{\mathit{k}}}$ on each element. For each element D^k a suitable set of local grid points ${{x}}_{\mathit{i}}^{\mathit{k}}\mathrm{\in }{\mathrm{D}}^{\mathit{k}}$ is chosen and the local continuous solution of ${\mathit{I}}_{\mathit{h}}^{\mathit{k}}$ is expressed using the interpolating Lagrange polynomials, such that ${\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{:}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{{x}}_{\mathit{i}}^{\mathit{k}}\mathrm{\right)}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)} \mathrm{forall}{x}\mathrm{\in }{\mathrm{D}}^{\mathit{k}}$(13)is a polynomial of order q on each element. Here the so called trial functions ${\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}$ are the interpolating two-dimensional Lagrange polynomials through the grid points ${{x}}_{\mathit{i}}^{\mathit{k}}$, i.e. ${\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{:}\mathrm{=}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}\mathit{x}\mathrm{\right)}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}\mathit{y}\mathrm{\right)}\mathrm{=}\underset{\begin{array}{c}\mathit{j}\mathrm{=}\mathrm{1}\\ \mathit{j}\ne \mathit{i}\end{array}}{\overset{{\mathit{N}}_{\mathit{p}}}{\mathrm{􏽙}}}\frac{\mathrm{(}\mathit{x}\mathrm{-}{\mathit{x}}_{\mathit{j}}\mathrm{)}}{\mathrm{(}{\mathit{x}}_{\mathit{i}}\mathrm{-}{\mathit{x}}_{\mathit{j}}\mathrm{)}}\hspace{0.17em}\underset{\begin{array}{c}\mathit{j}\mathrm{=}\mathrm{1}\\ \mathit{j}\ne \mathit{i}\end{array}}{\overset{{\mathit{N}}_{\mathit{p}}}{\mathrm{􏽙}}}\frac{\mathrm{(}\mathit{y}\mathrm{-}{\mathit{y}}_{\mathit{j}}\mathrm{)}}{\mathrm{(}{\mathit{y}}_{\mathit{i}}\mathrm{-}{\mathit{y}}_{\mathit{j}}\mathrm{)}}\mathrm{·}$(14)For a given order q the number of local grid points N_p is given by ${\mathit{N}}_{\mathit{q}}\mathrm{=}\frac{\mathrm{(}\mathit{q}\mathrm{+}\mathrm{1}\mathrm{)}\hspace{0.17em}\mathrm{(}\mathit{q}\mathrm{+}\mathrm{2}\mathrm{)}}{\mathrm{2}}$(15)for P_q elements and by ${\mathit{N}}_{\mathit{q}}\mathrm{=}\mathrm{(}\mathit{q}\mathrm{+}\mathrm{1}\mathrm{)}2$(16)for Q_q elements, respectively. In principle the polynomial order q and thus N_p can vary from element to element. However, for simplicity of notation we assume here that all elements have the same order q. Subsequently, the set of the local grid point indices within each element k is denoted by ${\mathrm{𝒩}}^{\mathit{k}}\mathrm{=}\mathrm{\left\{}\mathrm{1}\mathit{,}\mathit{...}\mathit{,}{\mathit{N}}_{\mathit{p}}\mathrm{\right\}}\mathit{.}$(17)The global piecewise continuous solution can then be written as the direct sum of all K solutions $\mathit{I}\mathrm{\left(}{x}\mathrm{\right)}\backsimeq {\mathit{I}}_{\mathit{h}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\underset{\mathit{k}\mathrm{=}\mathrm{1}}{\overset{\mathit{K}}{\mathrm{􏽍}}}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\mathit{.}$(18)The trial functions ${\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}$ introduced in Eq. (13)form a basis of the finite vector space V_h. This vector space is locally defined by ${\mathit{V}}_{\mathit{h}}^{\mathit{k}}\mathrm{=}\mathrm{spann}\mathrm{\left\{}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\right(}{\mathrm{D}}^{\mathit{k}}\mathrm{)}{\mathrm{\}}}_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}$ and globally given by ${\mathit{V}}_{\mathit{h}}\mathrm{=}{{}^{\mathrm{􏽌}}}_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\mathit{V}}_{\mathit{h}}^{\mathit{k}}$. Unlike the case of the continuous Galerkin finite element method (Eriksson et al. 1996), we note that no continuity of the trial functions (and therefore also of the solution I) is assumed or required across inter–element boundaries.

We define the residual of Eq. (10)by $\begin{array}{ccc}{\mathrm{ℛ}}_{\mathit{h}}\mathrm{\left(}{\mathit{I}}_{\mathit{h}}\mathrm{\right(}{x}\mathrm{\left)}\mathrm{\right)}\mathrm{=}\hspace{0.17em}& & \frac{\mathit{\partial }}{\mathit{\partial r}}\left(\mathit{\mu }\hspace{0.17em}{\mathit{I}}_{\mathit{h}}\right)\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}\right)\\ & & \mathrm{+}\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}\mathrm{+}\mathit{\chi ̂}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}\mathrm{-}{\mathit{\eta }}^{\mathrm{sp}}\mathrm{-}\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{I}}_{\mathit{h}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\end{array}$(19)with I_h ∈ V_h. In the context of finite element methods this residual is minimised with respect to all test functions ${\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{\in }{\mathit{C}}_{\mathrm{c}}^{\mathrm{\infty }}\mathrm{\left(}{\mathrm{D}}^{\mathit{k}}\mathrm{\right)}$ in each element. This gives the weak formulation of Eq. (10), which represents the physical principle of virtual work. Thus, it is required that ${\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}{\mathrm{ℛ}}_{\mathit{h}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}{\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}\mathrm{0} \mathrm{forall}{\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{\in }{\mathit{C}}_{\mathrm{c}}^{\mathrm{\infty }}\mathrm{\left(}{\mathrm{D}}^{\mathit{k}}\mathrm{\right)}\mathrm{and}\mathit{k}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,K}\mathit{.}$(20)For this study we consider only test functions, which are members of the same finite vector space as the trial functions, i.e. ${\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{\in }{\mathit{V}}_{\mathit{h}}^{\mathit{k}}\mathrm{\subset }{\mathit{C}}_{\mathrm{c}}^{\mathrm{\infty }}\mathrm{\left(}{\mathrm{D}}^{\mathit{k}}\mathrm{\right)}$. This corresponds to a Ritz-Galerkin method (Hesthaven & Warburton 2008). We note that in principle ${\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}$ may also be defined as a member of vector spaces other than ${\mathit{V}}_{\mathit{h}}^{\mathit{k}}$ (e.g. δ-distributions). The choice of different trial and test vector spaces leads to Petrov-Galerkin finite element schemes (Hesthaven & Warburton 2008) which are, however, not considered here.

For this particular choice of test functions we require for each element D^k${\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}{\mathrm{ℛ}}_{\mathit{h}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}{\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}\mathrm{0} \mathrm{forall}{\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{\in }{\mathit{V}}_{\mathit{h}}^{\mathit{k}}\mathrm{and}\mathit{k}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,K,}$(21)i.e. the residual is required to be orthogonal to all test functions ${\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{\in }{\mathit{V}}_{\mathit{h}}^{\mathit{k}}$. Since the test functions are given by a linear combination of the Lagrange polynomials, i.e. ${\mathit{\phi }}_{\mathit{h}}^{\mathit{k}}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}\mathit{\phi ̂}\begin{array}{c}\mathit{k}\\ \mathit{i}\end{array}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}$(22)with the expansion coefficients $\mathit{\phi ̂}\begin{array}{c}\mathit{k}\\ \mathit{h}\end{array}\mathrm{\in }R$, Eq. (21)results in the N_p equations ${\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}{\mathrm{ℛ}}_{\mathit{h}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}\mathrm{0} \mathrm{for}\mathit{i}\mathrm{=}\mathrm{1}\mathit{,}\mathit{...}\mathit{,N}p$(23)on each element D^k.

In summary, we consider in this work the following DG(q)-method: Seek I_h ∈ V_h such that $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\{\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}\mathit{\mu }\hspace{0.17em}{{\mathit{I}}_{\mathit{h}}^{\mathit{k}}}^{\mathrm{)}}\mathrm{+}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\left(\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\right)\mathrm{+}\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{+}\mathit{\chi ̂}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{-}{\mathit{\eta }}^{\mathrm{sp}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}\mathrm{-}\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}\mathrm{0}\begin{array}{c}\end{array}\mathrm{in}\mathrm{\Omega }\mathit{,}\\ \begin{array}{c}\end{array}\frac{\mathit{\partial }{\mathit{I}}_{\mathit{h}}}{\mathit{\partial r}}\mathrm{=}\mathrm{-}\frac{\mathrm{2}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}\mathrm{\mp }\mathit{\chi ̂}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}\mathrm{\pm }{\mathit{\eta }}^{\mathrm{sp}}\\ \begin{array}{c}\end{array}\mathrm{\pm }\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{I}}_{\mathit{h}}\mathrm{\left(}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{(}\mathrm{\pm }\mathrm{1}\mathit{,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\prime \begin{array}{c}\end{array}\mathrm{on}{\mathrm{\Gamma }}_{{\mathit{\mu }}^{\mathrm{\pm }}}\mathit{,}\\ \begin{array}{c}\end{array}{\mathit{I}}_{\mathit{h}}\mathrm{=}{\mathit{I}}_{\mathrm{in}}\begin{array}{c}\end{array}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{in}}^{\mathrm{+}}\mathit{,}\\ \begin{array}{c}\end{array}{\mathit{I}}_{\mathit{h}}\mathrm{=}{\mathit{I}}_{\mathrm{out}}\begin{array}{c}\end{array}\mathrm{on}{\mathrm{\Gamma }}_{\mathrm{out}}^{\mathrm{-}}\end{array}$(24)for i = 1,...,N_p and k = 1,...,K.

3.1. Elementwise calculations

Let us first consider the two terms $\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{+}\mathit{\chi ̂}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}$ under the integral in Eq. (24). Using the expansion of the solution ${\mathit{I}}_{\mathit{h}}^{\mathit{k}}$ within one element k (see Eq. (13)) yields $\begin{array}{ccc}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\left(\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\mathrm{+}\mathit{\chi ̂}\right)\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{{x}}_{\mathit{j}}^{\mathit{k}}\mathrm{\right)}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}& & \left(\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\mathrm{+}\mathit{\chi ̂}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{j}}^{\mathit{k}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\hspace{0.17em}\mathrm{d}{x}\\ & & \mathrm{forall}\mathit{i,j}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathit{.}\end{array}$(25)This defines the local symmetric mass matrix ${\mathrm{ℳ}}^{\mathit{k}}\mathrm{=}\mathrm{\left(}{\mathrm{ℳ}}_{\mathit{ij}}^{\mathit{k}}\mathrm{\right)}$ with ${\mathrm{ℳ}}_{\mathit{ij}}^{\mathit{k}}\mathrm{:}\mathrm{=}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\left(\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\mathrm{+}\mathit{\chi ̂}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{j}}^{\mathit{k}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{x} \mathrm{forall}\mathit{i,j}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathit{.}$(26)The term η^sp in Eq. (24)simply yields the local load vector ${\mathrm{ℬ}}^{\mathit{k}}\mathrm{=}\mathrm{\left(}{\mathrm{ℬ}}_{\mathit{i}}^{\mathit{k}}\mathrm{\right)}$, with ${\mathrm{ℬ}}_{\mathit{i}}^{\mathit{k}}\mathrm{:}\mathrm{=}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}{\mathit{\eta }}^{\mathrm{sp}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{x} \mathrm{forall}\mathit{i}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathit{.}$(27)Next we consider the differential with respect to r in Eq. (24). An integration by parts shifts the differential from the unknown quantity ${\mathit{I}}_{\mathit{h}}^{\mathit{k}}$ to the test function ${\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}$ introducing an additional integral over the boundary of the element D^k, i.e. $\begin{array}{ccc}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\frac{\mathit{\partial }}{\mathit{\partial r}}\mathrm{(}\mathit{\mu }\hspace{0.17em}{{\mathit{I}}_{\mathit{h}}^{\mathit{k}}}^{\mathrm{)}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{x}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\mathit{\mu }\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}}{\mathit{\partial r}}\mathrm{d}{x}\mathrm{+}{\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}& & \hat{{n}}\mathit{r}\mathrm{·}\mathit{\mu }\hspace{0.17em}{\mathit{I}}_{\mathit{h}}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{s}\\ & & \mathrm{forall}\mathit{i}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathit{,}\end{array}$(28)where $\hat{{n}}{}_{\mathit{r}}$ denotes the outward facing normal vector in r direction (see Fig. 3).

Because it is allowed to be discontinuous across element boundaries, the solution I_h is not uniquely defined at the boundary. Therefore, a numerical flux (f)^∗ is introduced which describes the flow of information from one element to another, yielding ${\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}\hat{{n}}\mathit{r}\mathrm{·}\mathit{\mu }\hspace{0.17em}{\mathit{I}}_{\mathit{h}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{s}\mathrm{=}{\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}\hat{{n}}\mathit{r}\mathrm{·}\mathrm{\left(}\mathit{f}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{\right(}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathit{,}{\mathit{I}}_{\mathit{h}}^{\mathit{n}}\mathrm{)}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{s}\mathit{.}$(29)The numerical flux depends in general on the solutions of the local element ${\mathit{I}}_{\mathit{h}}^{\mathit{k}}$ and the adjacent element ${\mathit{I}}_{\mathit{h}}^{\mathit{n}}$ and is directly related to the dynamics of the considered differential equation. For a particular choice of the flux a matrix ℱ^kn can be defined which contains the contributions due to the jump at the inter-element boundaries between the elements k and n. The exact contribution depends on the details of the numerical flux used, as discussed in Sect. 3.2.

Using again the expansion of Eq. (13)for the solution in the element k, the local stiffness matrix ${\mathrm{𝒮}}_{\mathit{r}}^{\mathit{k}}\mathrm{=}\mathrm{\left(}{\mathrm{𝒮}}_{\mathit{ij,r}}^{\mathit{k}}\mathrm{\right)}$ can be introduced as ${\mathrm{𝒮}}_{\mathit{ij,r}}^{\mathit{k}}\mathrm{:}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\mathit{\mu }\hspace{0.17em}{\mathit{\ell }}_{\mathit{j}}^{\mathit{k}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}}{\mathit{\partial r}}\hspace{0.17em}\mathrm{d}{x} \mathrm{forall}\mathit{i,j}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathit{.}$(30)For the differential with respect to μ in Eq. (24)integration by parts yields $\begin{array}{ccc}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\frac{\mathit{\partial }}{\mathit{\partial \mu }}\hspace{0.17em}& & \left(\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}\\ & & \mathrm{-}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\frac{\mathit{\partial }{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{d}{x}\mathrm{+}{\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}\hat{{n}}\mathit{\mu }\mathrm{·}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\hspace{0.17em}\mathrm{d}{s}\\ \mathrm{=}& & \mathrm{-}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{d}{x}\mathrm{+}{\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}\hat{{n}}\mathit{\mu }\mathrm{·}\mathrm{(}\mathit{f}{\mathrm{)}}^{\mathrm{\ast }}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\hspace{0.17em}\mathrm{d}{s}\end{array}$(31)for all and the respective outward facing normal vector in μ direction (see Fig. 3).

Fig. 3 Jump in the numerical solution along the normal $\hat{{n}}$ between P2 elements. The normals in the directions of r and μ are denoted by $\hat{{n}}{}_{\mathit{r}}$ and $\hat{{n}}{}_{\mathit{\mu }}$, respectively.

The corresponding stiffness matrix is given by ${\mathrm{𝒮}}_{\mathit{ij,\mu }}^{\mathit{k}}\mathrm{:}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{j}}^{\mathit{k}}\frac{\mathit{\partial }{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}}{\mathit{\partial \mu }}\hspace{0.17em}\mathrm{d}{x} \mathrm{forall}\mathit{i,j}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathit{.}$(32)The total local stiffness matrix ${\mathit{S}}^{\mathit{k}}\mathrm{=}\mathrm{\left(}{\mathrm{𝒮}}_{\mathit{ji}}^{\mathit{k}}\mathrm{\right)}$ for the element k can then be obtained by the sum of ${\mathit{S}}_{\mathit{ij,\mu }}^{\mathit{k}}$ and ${\mathit{S}}_{\mathit{ij,r}}^{\mathit{k}}$, i.e. ${\mathrm{𝒮}}_{\mathit{ij}}^{\mathit{k}}\mathrm{:}\mathrm{=}{\mathrm{𝒮}}_{\mathit{ij,r}}^{\mathit{k}}\mathrm{+}{\mathrm{𝒮}}_{\mathit{ij,\mu }}^{\mathit{k}}\mathit{.}$(33)

3.2. Numerical flux

The numerical flux considered here is the known as the upwind flux (e.g. Cockburn 2003; Hesthaven & Warburton 2008), given by $\mathrm{(}\mathit{f}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{=}\mathrm{(}\mathit{aI}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{=}\mathrm{\left\{}\mathrm{\right\{}\mathit{aI}\mathrm{\left\}}\mathrm{\right\}}\mathrm{+}\frac{\mathrm{1}}{\mathrm{2}}\left|\mathit{a}\right|\mathrm{˜}\mathit{I}\mathrm{¨}$(34)with $\mathrm{\left\{}\mathrm{\right\{}\mathit{I}\mathrm{\left\}}\mathrm{\right\}}\mathrm{=}\frac{{\mathit{I}}^{\mathrm{-}}\mathrm{+}{\mathit{I}}^{\mathrm{+}}}{\mathrm{2}}$(35)as the average of the solutions and $\mathrm{˜}\mathit{I}\mathrm{¨}\mathrm{=}\hat{{n}}\mathrm{-}{\mathit{I}}^{\mathrm{-}}\mathrm{+}\hat{{n}}\mathrm{+}\mathit{I}+$(36)as the jump of the solution along the normal. The situation is shown in Fig. 3, where the jump between two P₂ elements is illustrated. The parameter a is here either given by $\mathit{a}\mathrm{=}\mathit{\mu }$(37)for the differential with respect to r, or by $\mathit{a}\mathrm{=}\frac{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}{\mathit{r}}$(38)for the differential with respect to μ. This particular flux was chosen because for the radiative transfer problem considered here information propagates only in the direction of photon travel. The upwind flux for radiative transfer problems has already been considered by e.g. Li (2006). We note, however, that the choice of this flux is not unique. Other fluxes with different properties may be used (see Cockburn 2003, for a review on DG fluxes). The connection between the test and trial functions of the adjacent elements can be written as a matrix $\begin{array}{ccc}{\mathrm{ℱ}}_{\mathit{ij}}^{\mathit{kn}}\mathrm{:}\mathrm{=}{\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}\hat{{n}}\mathit{r}\mathrm{·}\mathrm{(}\mathit{f}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{(}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\right)}\mathit{,}{\mathit{I}}_{\mathit{h}}^{\mathit{n}}\mathrm{(}{\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}{\mathrm{)}}^{\mathrm{)}}& & \hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{x}\\ & & \end{array}$(39)Unlike the other matrices (such as ℳ^k), which are only defined over a single element, the matrix ℱ^kn connects a local element k with adjacent elements n. If the scattering integral considered below is not taken into account, this would be the only connection between different elements.

3.3. Scattering matrix

The contribution due to scattering in Eq. (24)is given by the integral ${\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\left(\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x} \mathrm{forall}\mathit{i}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathit{.}$(40)Using the expansions from Eqs. (13)and (18)results in $\begin{array}{ccc}\int {\mathrm{D}}^{\mathit{k}}\left(\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}\mathrm{=}& & \\ \sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{K}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{\mathit{r}}_{\mathit{j}}^{\mathit{n}}\mathit{,}{\mathit{\mu }}_{\mathit{j}}^{\mathit{n}}\mathrm{\right)}& & \\ \mathrm{\times }{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\left(\overline{\mathit{s}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\hspace{0.17em}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}& & \end{array}$(41)for all and . This defines a matrix containing the contributions due to scattering $\begin{array}{ccc}{\mathrm{𝒥}}_{\mathit{ij}}^{\mathit{kn}}\mathrm{=}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\left(\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}& & \\ \mathrm{forall}\mathit{i}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{k}}\mathrm{and}\mathit{j}\mathrm{\in }{\mathrm{𝒩}}^{\mathit{n}}\mathit{.}& & \end{array}$(42)We note that the scattering integral in Eq. (42)only contains known quantities, namely the scattering phase function p⁽⁰⁾ and the Lagrange polynomials ℓ_j. Thus, the integral can be evaluated exactly for any arbitrary complex phase function, independent of the mesh resolution. This is a direct advantage over other radiative transfer methods, for example the discrete ordinate methods (Thomas & Stamnes 2002) where the evaluation of the scattering integral is directly coupled to the resolution of the angular mesh (Chandrasekhar 1960).

3.4. Transformation matrix

To assemble a global matrix for the solution of the resulting system of equations, the matrixes and tensors with the local indices i and j must be mapped onto a matrix with global indices. This is done via a transformation matrix . Let g be a linear bijective mapping, which maps the set of local indices i within each element k onto a set of global coordinates $\mathrm{𝒩}\mathrm{=}\mathrm{\{}\mathrm{1}\mathit{,}\mathit{...}\mathit{,k}\mathrm{·}{\mathit{N}}_{\mathit{p}}\mathrm{\}}$(43)with $\mathit{g}\mathrm{:}{\mathrm{𝒩}}^{\mathit{k}}\mathrm{\to }\mathrm{𝒩}\mathit{,} \mathit{g}\mathrm{\left(}\mathit{i,k}\mathrm{\right)}\mathrm{=}\mathit{\alpha }\mathrm{\in }\mathrm{𝒩}\mathit{.}$(44)The linear mapping can be written as a transformation matrix . The matrix ${\mathrm{𝒯}}_{\mathit{\alpha i}}^{\mathit{k}}$ then maps, for example, the scattering matrix onto the global scattering matrix with the global indices α and β, which yields ${\mathrm{𝒥}}_{\mathit{\alpha \beta }}\mathrm{=}\sum _{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathit{K}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}\sum _{\mathit{j}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}{\mathrm{𝒯}}_{\mathit{\beta j}}^{\mathit{n}}\hspace{0.17em}{\mathrm{𝒯}}_{\mathit{\alpha i}}^{\mathit{k}}\hspace{0.17em}{\mathrm{𝒥}}_{\mathit{ij}}^{\mathit{kn}}\mathit{.}$(45)The global mass matrix ℳ, stiffness matrix , and jump matrix ℱ can be obtained via the same process. The total left-hand side matrix of the linear system of equations can be expressed by

${\mathrm{𝒜}}_{\mathit{\alpha \beta }}\mathrm{=}{\mathrm{𝒥}}_{\mathit{\alpha \beta }}\mathrm{+}{\mathrm{ℳ}}_{\mathit{\alpha \beta }}\mathrm{+}{\mathrm{𝒮}}_{\mathit{\alpha \beta }}\mathrm{+}{\mathrm{ℱ}}_{\mathit{\alpha \beta }}\hspace{0.17em}\mathit{.}$(46)The resulting linear system of equations is then $\mathrm{𝒜}\mathrm{·}{I}\mathrm{=}\mathrm{ℬ}$(47)with ${{I}}_{\mathit{\alpha }}\mathrm{:}\mathrm{=}{\sum }_{\mathit{k}\mathrm{=}\mathrm{1}}^{\mathit{K}}{\sum }_{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}{\mathrm{𝒯}}_{\mathit{\alpha i}}^{\mathit{k}}\hspace{0.17em}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\left(}{{x}}_{\mathit{i}}^{\mathit{k}}\mathrm{\right)}$ and ℬ the global load vector (cf. Eq. (27)). Since is a sparse matrix, the system can be efficiently solved by an iterative method (such as GMRES) and an appropriate pre-conditioner in many cases (Hackbusch 1994).

However, as mentioned by Kanschat & Meinköhn (2004) the matrix is sometimes badly conditioned, especially in cases with high opacity and single scattering albedos. The usual pre-conditioners (such as ILU, Hackbusch 1994) in combination with GMRES then fail to yield a converged solution. In this case a more adapted pre-conditioner is required to solve the system of Eq. (47). An example of such a pre-conditioner based on the Eddington approximation but limited to a plane-parallel situation is given by Kanschat & Meinköhn (2004). Nonetheless such a pre-conditioner should also be adaptable to spherically symmetric atmospheres which will be the subject of a forthcoming publication.

4. Test problems

In this section, we test our previously described numerical solution of the radiative transfer equation with a discontinuous Galerkin method. Due to the lack of exact solutions, we need to rely on test problems for comparison, where the resulting radiation field is known a priori, or to make a comparison with results of previously published calculations.

The first two tests evaluate the numerical stability and flux conservation in an empty atmosphere. In this case, the photons follow the characteristic paths of the transfer Eq. (1). These scenarios also feature discontinuities in their solutions, which will serve as a test for the numerical stability of our DG formulation.

For the third scenario, we refer to the radiative transfer calculations for an isotropically scattering, spherical atmosphere with a power-law extinction coefficient (Hummer & Rybicki 1971). We compare our results with their published values and also verify that the energy is also conserved in this case as required. We note that the intention of these calculations is to validate and verify the numerical approach and its capabilities to cope with discontinuities or highly scattering atmospheres, for example.

4.1. Computational details

We implemented the numerical algorithms described in the previous section by using the general finite element library GetFEM++¹. In principle, the element-wise calculations can also be performed in the sense of Alberty et al. (1999). Details of our numerical implementation are given in the Appendix.

A DG(1) method on a structured grid with rectangular Q₁ elements is used for the first two test problems. For both cases, the grid contains a total of 100 grid points in the radial direction and 80 grid points in the angular direction. The radius grid points are distributed equidistantly, while the angular mesh is constructed by using a Gaussian quadrature rule. The relatively high resolution is required to capture the discontinuous solutions properly. In principle, a much lower resolution could be used in these cases if an unstructured grid is constructed where the elements are directly aligned to capture the (previously) known discontinuities in the solution.

For the third test problem, we use a DG(2) method on a structured Q₂ grid with 25 radial points and 10 angular values. The radial nodes are placed in logarithmic equidistant steps, while for the angular grid a Gaussian quadrature rule is again employed.

The resulting linear system of equations is solved by an iterative GMRES solver with a restart value of 10 in each case. The iteration matrix is preconditioned by an incomplete LU factorisation.

4.2. Empty atmospheres with an isotropically emitting inner surface

Fig. 4 Radiation field I(μ,r) as a function of radius r and angular variable μ for an empty, extended atmosphere with an isotropically emitting inner surface.

As a very first test we study the following simple radiative transfer problem $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{r}}_{\mathrm{in}}\mathrm{=}\mathrm{1.0}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{r}}_{\mathrm{out}}\mathrm{=}\mathrm{3.0}\mathit{,}\\ \begin{array}{c}\end{array}\mathit{\kappa }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\eta }}^{\mathrm{sp}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\eta }}^{\mathrm{ind}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\overline{\mathit{s}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{0}\end{array}$(48)with boundary conditions $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{I}}_{\mathrm{in}}\mathrm{=}\mathrm{4}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{I}}_{\mathrm{out}}\mathrm{=}\mathrm{0}\mathit{,}\end{array}$(49)which represents an empty, extended atmosphere with an isotropically emitting inner surface.

Fig. 5 Radiation field I(μ,r) as a function of the angular variable μ and three chosen radii r. The vertical dotted line denotes μc, the limiting value of μ under which the inner emitting surface is seen at the outer boundary r = 3. The discontinuity is accurately captured (i.e. the solution shows no oscillations).

The resulting radiation field I(μ,r) is shown in Fig. 4. Additionally, several cuts through the parameter space are given in Fig. 5, where the intensity I is shown as a function of μ for three different values of r. We note that in this radiative transfer problem a real discontinuity is present in the solution which separates the region of constant intensity originating from the inner surface and the regions where the intensity is zero. The characteristic separating these two regions is given by ${\mathit{\mu }}_{\mathrm{c}}\mathrm{=}\sqrt{\mathrm{1}\mathrm{-}\frac{{\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}}{{\mathit{r}}^{\mathrm{2}}}}\mathrm{·}$(50)The results in Figs. 4 and 5 clearly show the spherical dilution of the outward directed intensity. With increasing radius r the angle μ_c, under which the central emitting surface can be seen, decreases strongly. This is illustrated in greater detail in Fig. 5, where the intensity I is shown as a function of μ and for three different values of r. The vertical dotted line marks the limiting value of μ_c(r_out) under which the emitting inner surface is seen at the outer boundary. As suggested by Fig. 5, the DG method presented here accurately reproduces this spherical dilution of the radiation field without introducing numerical oscillations.

Consequently, the DG method is capable of capturing this discontinuity quite well without introducing numerical oscillations which would occur in case of a continuous Galerkin method. In practice, the discontinuity is smeared out over one element. In order to limit the errors introduced by this out-smearing it is possible to use either high mesh resolution in general or adaptive mesh refinement along the discontinuous boundaries of the intensity given by Eq. (50).

Because there are no interactions of the atmosphere with the radiation field in this test example, the radiative transfer problem should be automatically energy conserving. For a spherically symmetric atmosphere this means that the flux (r²H(r)) is constant throughout the atmosphere, where H is the first angular moment (Eddington flux) of the intensity $\mathit{H}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathit{I}\mathrm{\left(}\mathit{r,\mu }\mathrm{\right)}\hspace{0.17em}\mathit{\mu }\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathit{.}$(51)In this first test calculation, the flux is conserved better than 0.01% everywhere in the computational domain, indicating that the DG method is stable and accurate for such a problem.

4.3. Empty atmospheres illuminated by an exterior light source

This test case differs from the previous one only by choosing different boundary conditions $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{I}}_{\mathrm{in}}\mathrm{=}\mathrm{0}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{I}}_{\mathrm{out}}\mathrm{=}\mathrm{4}\mathit{\mu }\mathit{.}\end{array}$(52)This problem then represents an empty, extended atmosphere illuminated by an exterior light source. The particular form of the outer boundary condition has been chosen to highlight the photons’ characteristics, i.e. the geometric pathways of photons travelling through the atmosphere. The resulting radiation field is shown in Figs. 6 and 7.

Fig. 6 Radiation field I(μ,r) as a function of radius r and angular variable μ for an empty, extended atmosphere illuminated by an exterior light source. The curved black lines illustrate a number of selected characteristics of this radiative transfer problem, i.e. the pathways of photons through the atmosphere. The vertical white line separates ingoing and outgoing directions.

Fig. 7 Radiation field I(μ,r) as a function of the angular variable μ and three chosen radii r. The discontinuity is accurately captured (i.e. there are no oscillations in the solution).

As can be clearly be seen in Fig. 6, only a small fraction of the incident radiation eventually reaches the inner surface. The remaining radiation is transmitted through the atmosphere and leaves the computational domain at the outer boundary with the same intensity and under the same angle as the incident radiation field. The curved black lines in Fig. 6 illustrate some of the characteristics of this radiative transfer problem, i.e. the pathways of photons through the atmosphere. This behaviour is very different to problems with plane-parallel geometry where the incident light under any angle would eventually reach the surface.

Like the previous test case, this radiative transfer problem also has a discontinuous solution in a part of its computational domain. This discontinuity is also captured accurately by the DG formulation introduced in this work (see Fig. 7). Again, the flux is conserved better than 0.01% everywhere in the computational domain.

4.4. Isotropically scattering atmosphere

To prove that our DG radiative transfer method can handle non-local scattering problems, we compare our results with calculations published in Hummer & Rybicki (1971) resembling a purely scattering, spherical atmosphere with a power-law opacity in this section. In contrast to our work, Hummer & Rybicki (1971) used an iterative moment method to obtain the solution of the radiative transfer equation, although owing to the limitation of this method, only isotropic scattering was taken into account. The radiative transfer problem considered here for comparison is therefore specified by $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{r}}_{\mathrm{in}}\mathrm{=}\mathrm{0.01}\mathit{,}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{r}}_{\mathrm{out}}\mathrm{=}\mathrm{0.1}\mathit{,}\\ \begin{array}{c}\end{array}\mathit{\kappa }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\eta }}^{\mathrm{sp}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{\eta }}^{\mathrm{ind}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}\\ \begin{array}{c}\end{array}\overline{\mathit{s}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}{\mathit{r}}^{\mathrm{-}\mathrm{3}\mathit{/}\mathrm{2}}\mathit{.}\end{array}$(53)Following Hummer & Rybicki (1971), scattering is assumed to be isotropic, i.e. the scattering phase function is given by ${\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}\mathrm{1.}$(54)Boundary conditions are chosen to resemble those introduced in Hummer & Rybicki (1971). The moment method employed by Hummer & Rybicki (1971) allowed for a direct prescription of the radiation flux at the inner boundary, which was set to ${\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\mathit{H}\mathrm{\left(}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}\mathrm{=}\mathrm{1}$. At the outer boundary, no incident radiation was assumed. To replicate these constant flux boundary conditions, we employed the following conditions on the specific intensity at the domain boundaries $\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}{\mathit{I}}_{\mathrm{in}}\mathrm{=}\mathrm{4}\hspace{0.17em}{\mathit{r}}_{\mathrm{in}}^{-2}\mathrm{-}\mathrm{2}{\mathrm{\int }}_{-1}^{\mathrm{0}}\mathit{\mu }\hspace{0.17em}\mathit{I}\mathrm{\left(}\mathit{\mu ,}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array}{\mathit{I}}_{\mathrm{out}}\mathrm{=}\mathrm{0.}\end{array}$(55)This inner boundary condition is obtained from the requirement that $\begin{array}{ccc}{\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\mathit{H}\mathrm{\left(}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}& \mathrm{=}& {\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}\mathit{\mu }\hspace{0.17em}\mathit{I}\mathrm{\left(}\mathit{\mu ,}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\\ & \mathrm{=}& {\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{0}}\mathit{\mu }\hspace{0.17em}\mathit{I}\mathrm{\left(}\mathit{\mu ,}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathrm{+}{\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{+}\mathrm{1}}\mathit{\mu }\hspace{0.17em}\mathit{I}\mathrm{\left(}\mathit{\mu ,}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\\ & \mathrm{=}& {\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\frac{\mathrm{1}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{0}}\mathit{\mu }\hspace{0.17em}\mathit{I}\mathrm{\left(}\mathit{\mu ,}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}\hspace{0.17em}\mathrm{d}\mathit{\mu }\mathrm{+}{\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\hspace{0.17em}\frac{\mathrm{1}}{\mathrm{4}}\hspace{0.17em}\mathit{I}\mathrm{(}\mathrm{+}\mathit{\mu ,}{\mathit{r}}_{\mathrm{in}}\mathrm{)}\\ & \mathrm{=}& \mathrm{1}\mathit{,}\end{array}$(56)assuming an isotropic distribution of I( + μ,r_in).

Fig. 8 Radiation field I(μ,r) as a function of radius r and angular variable μ. The vertical black line separates ingoing and outgoing directions.

Because of the purely scattering atmosphere, the radiation flux r²H(r) needs to be conserved throughout the atmosphere. The resulting values of r²H(r) differ from the assumed constant flux at the inner boundary ${\mathit{r}}_{\mathrm{in}}^{\mathrm{2}}\mathit{H}\mathrm{\left(}{\mathit{r}}_{\mathrm{in}}\mathrm{\right)}\mathrm{=}\mathrm{1}$ by less than 0.5% everywhere in the computational domain. In the inner part of the atmosphere, the highest deviations are below 0.001% indicating again the capability of the DG method.

The resulting radiation field I(μ,r) is shown in Fig. 8. In contrast to the previous two test problems, no discontinuity occurs in this case. Owing to the many multiple scattering events, the radiation field exhibits a smooth behaviour in angle and radius. At the inner boundary, the radiation field is isotropic, resulting in an Eddington factor of 0.3 as expected.

Mean intensity.

Figure 9 shows the resulting values of the mean intensity r²J(r) in comparison to the corresponding tabulated values taken from Table II in Hummer & Rybicki (1971). Agreement with the results of Hummer & Rybicki (1971) are excellent in all parts of the isotropically scattering atmosphere. It should be noted that the grid resolution used in the present calculation is much lower than that employed by Hummer & Rybicki (1971) (200 r-points, 100 angles). We use a factor of eight fewer radius grid points and a factor of ten fewer angular values. The radiation field as a function of radius and angle is rather smooth and so does not require a high resolution – especially since the intensity here is assumed to be a second-order polynomial within each element.

Fig. 9 Mean intensity r2J as a function of the radius for an isotropically scattering atmosphere. The dots mark the corresponding results published by Hummer & Rybicki (1971).

Limb-darkening.

The limb-darkening curves are shown in Fig. 10 as a function of the impact parameter , given by $\mathrm{𝒫}\mathrm{=}{\mathit{r}}_{\mathrm{out}}\sqrt{\mathrm{1}\mathrm{-}{\mathit{\mu }}^{\mathrm{2}}}\mathit{.}$(57)The curve shows an excellent agreement with the results of Hummer & Rybicki (1971; their Fig. 6). The calculated disk centre intensity I(0) is 816.4, which deviates by less than 0.5% from the value of 820 given in Hummer & Rybicki (1971). We note, however, that their solution was stated to have a flux conservation error of a few percent at the outer boundary, while our flux is conserved better than 0.5% everywhere, which can explain the slight deviations in the disk centre intensity.

Fig. 10 Limb-darkening curve I(p) /I(0) as a function of the impact parameter . The dots mark the corresponding results from Hummer & Rybicki (1971), which have been directly read off their Fig. 6.

5. Numerical efficiency

In this section we study the development of the computational time as a function of the total number of degrees of freedom N by using the previous, third test example for the DG(2) method (isotropic scattering atmosphere). We test separately the dependence on the number of angular and on radius points. Successively increasing the resolution in either grid direction, we track the development of the computational effort. For a grid composed of Q_q elements with N_q local degrees of freedom N is given by $\mathit{N}\mathrm{=}\mathrm{(}{\mathit{N}}_{\mathit{r}}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{·}\mathrm{(}{\mathit{N}}_{\mathit{\mu }}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{·}{\mathit{N}}_{\mathit{q}}\mathit{,}$(58)where N_r and N_μ are the number of radial and angular grid points. According to ${\mathit{N}}_{\mathit{q}}\mathrm{=}\mathrm{(}\mathit{q}\mathrm{+}\mathrm{1}{\mathrm{)}}^{\mathrm{2}}\mathit{,}$(59)the number of local degrees of freedom for the Q₂ elements used here is nine (see Fig. 2). We note that the computational time includes every single step to solve the problem, from the construction of the grid and finding the corresponding elements throughout the grid that need to be connected by the scattering integral to calculating the element-wise contributions to the matrix and solving the linear system of equations. In reality, many of these tasks need to be done only once – such as the element connections – and their results can be re-used at different iteration steps or for other frequencies.

The corresponding results are shown in Fig. 11, which also depicts the corresponding linear fits through the obtained data points to determine the order of dependence on the grid point number. The results in Fig. 11 suggest that the computing time increases linearly with the number of radius points. For a smaller number of angles, the computational time increases with N^1.5. For numbers of N larger than about 10 000, this dependency changes to N^2.5 (see Fig. 11). For comparison, the iterative method using variable Eddington factors as employed by Hummer & Rybicki (1971) typically scales cubic in the number of radius grid points for one iteration step Mihalas (1978), if for each radius point a corresponding tangent ray (equal to two angular directions) is used. Several iterations, however, are needed to obtain a converged solution, which means that the overall computational cost is larger than this estimate.

It should be noted that with increasing numbers of degrees of freedom, most of the time is spent solving the linear system of equations. This is, in fact, the main reason for the increase in the computational time for the angular points. The scattering integral (Eq. (42)) is responsible for connecting elements in the final matrix which are, however, not direct neighbours in the grid. Depending on the actual implementation of the element numbering scheme (and, thus their position in ), the contributions to the integral Eq. (42)by those elements can be scattered throughout the matrix . This makes an efficient solution of the linear system of equations quite hard for an iterative solver. A considerable decrease in computing time could, therefore, be obtained when the element numbering scheme would take the elements’ connections due to the scattering integral into account, such that the contributions of these elements would be located closer to the main diagonal of the matrix (see e.g. Cuthill & McKee 1969). For a purely scattering atmosphere, we were able to obtain converged solutions up to optical depths of around 100. Beyond that, the matrix becomes too badly conditioned to be solved by a standard GMRES scheme with an ILU preconditioner. As mentioned in Sect. 3.4, a more sophisticated preconditioner needs to be developed for such cases.

Fig. 11 Development of the normalised computational time as a function of the number of degrees of freedom. Both variables are depicted in logarithmic scale. The black dots mark the results for increasing the number of radius grid points from 40 to 500, while the number of angles is kept constant at 10. The red crosses indicate the results for increasing the angular resolution from 10 to 82, using a constant radial grid point number of 40. The corresponding solid lines depict the linear fits through the respective data points. For the angular nodes, two different fits are used.

6. Outlook

The implementation of the DG-FEM presented in this work can be greatly extended in the future in view of its numerical efficiency. One improvement which can be addressed is adaptive grid refinement in combination with error estimation. To this end, the residual in Eq. (23)provides the opportunity of obtaining an a posteriori error estimate. The pointwise error of the calculated solution can be estimated by expressing this equation in a suitable norm. In regions where the error is above a certain threshold, a local grid refinement or a local increase of the polynomial order can then be used to increase the numerical accuracy. Using this adaptive grid refinement ensures that the exact solution is calculated accurately on the entire grid. Adaptive mesh refinement can become important when the solution features a non-monotonous intensity distribution. This can occur within circumstellar dust shells, for example, where the dust forming region somewhere in the atmosphere creates a local maximum in the angular intensity distribution (Gail & Sedlmayr 2014). Such a non-monotonous behaviour can easily be traced by local grid refinement without increasing the resolution in other parts of the grid.

Another option for improving the numerical efficiency of the presented DG-FEM is to parallelise the numerical work. In principle, since the DG-FEM formulation does not enforce continuity across element boundaries, all element-wise calculations presented in Sect. 3.1 can be computed independently from each other. Thus, such independent computations can be parallelised very easily by using MPI or OpenMP schemes, for example. Furthermore, the calculations within each element – unless pushed to very high orders – are always a comparatively small numerical problem. Such small problems are ideal for running on GPU clusters, i.e. graphic card processors with thousands of cores, which could increase the computational speed by at least one order of magnitude.

Finally, the solution of the final linear system of equations can also be parallised. Special libraries, such as PETSc, are available to distribute the computational task across several CPUs or GPUs. This would be especially helpful in cases where most of the time is spent by the iterative solver.

7. Summary

In this work we presented a numerical method to directly solve the radiative transfer equation in spherical symmetry. We used a discontinuous Galerkin finite element method to discretise the transport equation as a two-dimensional boundary value problem. Arbitrary complex scattering phase function can, in particular, be integrated exactly using this approach. The discretisation yields a sparse system of equations which can be solved by standard iterative methods with suitable preconditioners.

The applicability of the DG-FEM is verified on a set of different test problems. The discontinuous Galerkin method is able to directly describe the spherical dilution of the radiation field in an extended atmosphere. We show that the DG-FEM is able to accurately capture discontinuities which might appear in the solution of the transfer equation for special test problems. We also compared our method with previously published radiative transfer calculations for a spherical, isotropically scattering atmosphere. The agreement between the different approaches was found to be excellent. Consequently, the discontinuous Galerkin method is perfectly suitable for treating accurately the radiative transfer in a spherical atmosphere even with, for example, discontinuities in the radiation field or complex scattering behaviour, which might cause severe numerical problems for other radiative transfer solution methods.

Acknowledgments

D.K. gratefully acknowledges the support of the Center for Space and Habitability of the University of Bern and the MERAC Foundation for partial financial assistance. This work has been carried out within the frame of the National Centre for Competence in Research PlanetS supported by the Swiss National Science Foundation. D.K. acknowledges the financial support of the SNSF. The authors thank the anonymous referee for the very constructive and helpful comments that greatly improved the manuscript.

References

Appendix A: Implementation details

We briefly give details of our implementation of the DG-FEM in a numerical code as one example. As previously mentioned, we use the general finite element library GetFEM++. This library provides a large set of high-level functions and tools which can assist in easily setting up a finite element method. A complete description of the library can be found in its user documentation². In this appendix we only include the basic set of necessary steps. We refer the reader to the GetFem++ documentation for much more detailed descriptions and explanations on the contents and implementations of the library.

If, alternatively, one is interested in building a finite element from scratch, we refer to e.g. Hesthaven & Warburton (2008), Eriksson et al. (1996), or Alberty et al. (1999) for more details.

Appendix A.1: Grid and element generation

The first step in setting up a finite element method is the generation of the grid and the definition of the elements. Within GetFEM++, the mesh generation can be done via different approaches. One possibility is to construct a grid outside of GetFEM++ (for example by using the mesh generator Gmsh) and to import it. It is also possible to use GetFEM++ itself for grid generation; however, the built-in mesh generation tools cover only very simple grids. A third possibility is to construct a special grid by defining all vertices of the mesh elements. We use this last option since, for example, we want to distribute the angular vertices according to a Gaussian quadrature rule which is not covered by the mesh tools provided by GetFEM++.

A single point is added to a mesh cell with

pts[i] = bgeot::base_node(radius_point,angular_point); ,

followed by

mesh.add_parallelepiped_by_points(2, pts.begin());

to create, e.g., a square mesh cell.

For a given mesh, GetFEM++ will then automatically create the elements. It offers a wide choice of different element definitions (see the user documentation for a complete overview). For our DG-FEM code, we use Lagrangian elements (see Fig. 2), which are created via

mesh_fem.set_finite_element(getfem::fem_descriptor("FEM_QK_DISCONTINUOUS(2,2)"));

for a second-order method on a two-dimensional grid. This has to be coupled with a suitable quadrature method, e.g.

mesh_integration.set_integration_method(getfem::int_method_descriptor ("IM_GAUSS_PARALLELEPIPED(2,4)"));

GetFem++ also offers other quadrature methods which are extensively described in the library’s documentation.

Given a mesh and the definition of the corresponding elements of a chosen polynomial degree, GetFEM++ automatically produces the transformation matrix ${\mathrm{𝒯}}_{\mathit{\alpha i}}^{\mathit{k}}$ to map the local coordinates i within the kth element to global coordinates α. GetFEM++ also provides tools to perform calculations of base functions and their derivatives, and to do computations on the faces of elements and providing the corresponding normal vectors at the elements’ faces.

Appendix A.2: Mass matrix

Within GetFEM++, all computations are usually defined on a reference element. The library provides high-level methods for the assembly of the mass and stiffness matrices. For example, the computation of the mass matrix on a single element (cf. Eq. (26)) ${\mathrm{ℳ}}_{\mathit{ij}}^{\mathit{k}}\mathrm{:}\mathrm{=}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\left(\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\mathrm{+}\mathit{\chi ̂}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{j}}^{\mathit{k}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{x}$(A.1)is simply defined by

"a=data(#1);" "M(#1,#1) += a(i).comp(Base(#1).Base(#1).Base(#1))(i,:,:)",

where Base(#1) are the base functions ℓ^k and a contains the function $\mathit{a}\mathrm{=}\frac{\mathrm{2}\hspace{0.17em}\mathit{\mu }}{\mathit{r}}\mathrm{+}\mathit{\chi ̂}\mathit{b,}$(A.2)which, following Eq. (13), is expanded into $\mathit{a}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathit{p}}}\mathit{a}\mathrm{\left(}{{x}}_{\mathit{i}}^{\mathit{k}}\mathrm{\right)}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\mathit{.}$(A.3)This summation is denoted i in the definition (i,:,:). The two remaining free indices (denoted by :) are the corresponding matrix indices.

The global mass matrix is defined according to

mass_matrix.set("a=data(#1);" "M(#1,#1) += a(i).comp(Base(#1).Base(#1).Base(#1))(i,:,:)");mass_matrix.assembly(); ,

where the local contributions of each element are calculated and then automatically mapped to the real elements by means of a Jacobian matrix.

Appendix A.3: Stiffness matrix

By analogy, the local stiffness matrix (cf. Eq. (30)) ${\mathrm{𝒮}}_{\mathit{ij,r}}^{\mathit{k}}\mathrm{:}\mathrm{=}\mathrm{-}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\mathit{\mu }\hspace{0.17em}{\mathit{\ell }}_{\mathit{j}}^{\mathit{k}}\hspace{0.17em}\frac{\mathit{\partial }{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}}{\mathit{\partial r}}\hspace{0.17em}\mathrm{d}{x}$(A.4)is defined by

"M(#1,#1) -= comp(Base(#2).Base(#1).Grad(#1))(1,:,:,1);".

The parameter Base(#2) describes the μ factor in Eq. (A.4), which in the context of GetFem++ is expressed by a global function. This function returns the μ value for any given point x. The derivative of the base function with respect to r is denoted Grad(#1). The 1 in the last entry in the definition (1,:,:,1) signifies that the derivative has to be performed with respect to the first dimension of the two-dimensional grid, which in our case is the radius coordinate. The global stiffness matrix is then assembled with

stiffness_matrix.set("M(#1,#1) -= comp(Base(#2).Base(#1).Grad(#1))(1,:,:,1);");stiffness_matrix.assembly(); . Of course, the second stiffness matrix (Eq. (32)) is calculated analogously.

Appendix A.4: Numerical flux

The numerical flux (Eq. (39)) is evaluated by integrating along the faces of each element, taking into account the contributions of adjacent elements: ${\mathrm{ℱ}}_{\mathit{ij}}^{\mathit{kn}}\mathrm{:}\mathrm{=}{\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}\hat{{n}}\mathit{r}\mathrm{·}\mathrm{(}\mathit{f}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{\left(}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\right(}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left)}\mathit{,}{\mathit{I}}_{\mathit{h}}^{\mathit{n}}\mathrm{\right(}{\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}{\mathrm{\left)}}^{\mathrm{\right)}}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{x}\mathit{.}$(A.5)Since the high-level methods from the previous section only operate on one local element, it is necessary to use the special getfem::compute_on_inter_element class which also includes information on the properties of elements sharing a common face. The class method compute_on_gauss_point has to be carefully adapted to calculate the contribution to the integral of Eq. (39)at each of the integral’s quadrature nodes. Within this special object class, GetFEM++ provides all necessary information to evaluate the integral, such as the normal vectors, quadrature nodes, and weights along the face, as well as the entries of the Jacobian matrix to map the computations on the reference element back to the real one. More details on the implementation of this class can be found in the official documentation and the corresponding code examples published together with the source code of the library.

Appendix A.5: Scattering matrix

The scattering matrix is given by (cf. Eq. (42)) ${\mathrm{𝒥}}_{\mathit{ij}}^{\mathit{kn}}\mathrm{=}{\mathrm{\int }}_{{\mathrm{D}}^{\mathit{k}}}\left(\frac{\overline{\mathit{s}}}{\mathrm{2}}{\mathrm{\int }}_{-1}^{\mathrm{+}\mathrm{1}}{\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}{\mathit{p}}^{\mathrm{\left(}\mathrm{0}\mathrm{\right)}}\mathrm{\left(}\mathit{r,\mu ,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{\mathit{\mu }}^{\mathrm{\prime }}\right)\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left(}{x}\mathrm{\right)}\hspace{0.17em}\mathrm{d}{x}\mathit{.}$(A.6)It contains two integrals, the outer one for the integration over a local element and the inner one which integrates the scattering phase function over multiple different elements. For this more complicated term GetFEM++ also provides no direct high-level method.

We use a Gaussian quadrature rule for the discretisation of the inner integral, with an order high enough to ensure an exact integration of the phase function. At each of the Gaussian quadrature nodes from this scattering integral μ′, the value of the base function ${\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}$ is needed; this point usually lies outside of the considered, local element k. Unless a degree of freedom is directly located at the coordinates (r,μ′), the base function ${\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}\mathrm{\left(}\mathit{r,}{\mathit{\mu }}^{\mathrm{\prime }}\mathrm{\right)}$ contains contributions from several different degrees of freedom within the element n.

GetFEM++ provides interpolation routines to find the corresponding values of ${\mathit{\ell }}_{\mathit{j}}^{\mathit{n}}$ for a given radius r and angle μ′, as well as the set of degrees of freedom which contribute to this specific base function. Given the mapping from the local coordinates of the ith degree of freedom in the kth element to the global coordinate α and the corresponding mappings of all degrees of freedom connected by the scattering integral with global indices β, it is possible to directly add their contributions to the global matrix .

Appendix A.6: Boundary conditions

The boundary conditions from Eq. (24)can be implemented in two different ways. They can either be introduced in analogy to the integral of the numerical jump from Eq. (39), where the intensity from the adjacent element is replaced by the boundary condition. For example, the outer boundary condition I_out can then be written as ${\mathrm{ℱ}}_{\mathit{ij}}^{\mathit{k}}\mathrm{=}{\mathrm{\int }}_{\mathit{\partial }{\mathrm{D}}^{\mathit{k}}}\hat{{n}}\mathit{r}\mathrm{·}\mathrm{(}\mathit{f}{\mathrm{)}}^{\mathrm{\ast }}\mathrm{\left(}{\mathit{I}}_{\mathit{h}}^{\mathit{k}}\mathrm{\right(}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{\left)}\mathit{,}{\mathit{I}}_{\mathrm{out}}\mathrm{\right)}\hspace{0.17em}{\mathit{\ell }}_{\mathit{i}}^{\mathit{k}}\mathrm{d}{x} {x}\mathrm{\in }{\mathrm{\Gamma }}_{\mathrm{out}}^{\mathrm{-}}\mathit{.}$(A.7)This, however, means that the boundary condition is only satisfied in an average sense because it is introduced as a variational formulation.

The second possibility is to directly manipulate the global matrix and the load vector ℬ. For any degree of freedom α located at an element’s face where a boundary condition needs to be applied, all entries are set to zero, while is set to unity. The αth component of the load vector is then set to the respective boundary condition. This ensures that the boundary conditions are always met exactly. For our implementation, we chose this second variant.

Appendix A.7: Solving the system of equations

GetFEM++ itself provides some basic iterative, numerical methods to solve the system of Eq. (47)optimised for sparse matrices, including a small set of preconditioners. In this study, we use the GMRES method with an ILU preconditioner.

However, if required, the sparse matrix and the load vector ℬ can also be forwarded to an external solver, such as PETSc or MUMPS, a parallel sparse direct solver library, for example.

All Tables

All Figures

	Fig. 1 Computational domain Ω. The boundaries of the domain between the inner radius rin and outer radius rout, as well as μ between + 1 and −1, are denoted by the respective Γ.
In the text

	Fig. 2 Lagrangian reference elements for dimension n = 2, and polynomial degrees q = 0,1,2,3. The degrees of freedom (see Table 1) are marked with •. Top: Pq elements, bottom: Qq elements.
In the text

	Fig. 3 Jump in the numerical solution along the normal $\hat{{n}}$ between P2 elements. The normals in the directions of r and μ are denoted by $\hat{{n}}{}_{\mathit{r}}$ and $\hat{{n}}{}_{\mathit{\mu }}$, respectively.
In the text

	Fig. 4 Radiation field I(μ,r) as a function of radius r and angular variable μ for an empty, extended atmosphere with an isotropically emitting inner surface.
In the text

	Fig. 5 Radiation field I(μ,r) as a function of the angular variable μ and three chosen radii r. The vertical dotted line denotes μc, the limiting value of μ under which the inner emitting surface is seen at the outer boundary r = 3. The discontinuity is accurately captured (i.e. the solution shows no oscillations).
In the text

	Fig. 6 Radiation field I(μ,r) as a function of radius r and angular variable μ for an empty, extended atmosphere illuminated by an exterior light source. The curved black lines illustrate a number of selected characteristics of this radiative transfer problem, i.e. the pathways of photons through the atmosphere. The vertical white line separates ingoing and outgoing directions.
In the text

	Fig. 7 Radiation field I(μ,r) as a function of the angular variable μ and three chosen radii r. The discontinuity is accurately captured (i.e. there are no oscillations in the solution).
In the text

	Fig. 8 Radiation field I(μ,r) as a function of radius r and angular variable μ. The vertical black line separates ingoing and outgoing directions.
In the text

	Fig. 9 Mean intensity r2J as a function of the radius for an isotropically scattering atmosphere. The dots mark the corresponding results published by Hummer & Rybicki (1971).
In the text

	Fig. 10 Limb-darkening curve I(p) /I(0) as a function of the impact parameter . The dots mark the corresponding results from Hummer & Rybicki (1971), which have been directly read off their Fig. 6.
In the text

Fig. 11 Development of the normalised computational time as a function of the number of degrees of freedom. Both variables are depicted in logarithmic scale. The black dots mark the results for increasing the number of radius grid points from 40 to 500, while the number of angles is kept constant at 10. The red crosses indicate the results for increasing the angular resolution from 10 to 82, using a constant radial grid point number of 40. The corresponding solid lines depict the linear fits through the respective data points. For the angular nodes, two different fits are used.

In the text