Issue |
A&A
Volume 690, October 2024
|
|
---|---|---|
Article Number | A192 | |
Number of page(s) | 14 | |
Section | Numerical methods and codes | |
DOI | https://doi.org/10.1051/0004-6361/202244322 | |
Published online | 09 October 2024 |
Discontinuous Galerkin finite element method for the continuum radiative transfer problem inside axis-symmetric circumstellar envelopes
Université de la Côte d’Azur, Observatoire de la Côte d’Azur, CNRS,
Laboratoire Lagrange, Bd de l’Observatoire, CS 34229,
06304
Nice cedex 4,
France
★ Corresponding author; e-mail: jeremy.perdigon@oca.eu
Received:
22
June
2022
Accepted:
16
August
2024
Context. The study of the continuum radiative transfer problem inside circumstellar envelopes is both a theoretical and numerical challenge, especially in the frequency-dependent and multi-dimensional case. While approximate methods are easier to handle numerically, they often fail to accurately describe the radiation field inside complex geometries. For these cases, it is necessary to directly solve the radiative transfer equation numerically.
Aims. We investigate the accuracy of the discontinuous Galerkin finite element method (DGFEM hereafter) applied to the frequency-dependent two-dimensional radiative transfer problem, and coupled with the radiative equilibrium equation. We next used this method in the context of axis-symmetric circumstellar envelopes.
Methods. The DGFEM is a variant of finite element methods. It employs discontinuous elements and flux integrals along their boundaries, ensuring local flux conservation. However, as opposed to the classical finite element methods, the solution is discontinuous across element edges. We implemented this approach in a code and tested its accuracy by comparing our results with the benchmarks from the literature.
Results. For all the tested cases, the temperatures profiles agree within one percent. Additionally, the emerging spectral energy distributions (SEDs) and images, obtained by ray-tracing techniques from the DGFEM emissivity, agree on average within 5% and 10%, respectively.
Conclusions. We show that the DGFEM can accurately describe the continuum radiative transfer problem inside axis-symmetric circumstellar envelopes. Consecutively the emerging SEDs and images are also well reproduced. The DGFEM provides an alternative method (other than Monte-Carlo methods for instance) for solving the radiative transfer equation, and it could be used in cases that are more difficult to handle with the other methods.
Key words: radiative transfer / methods: numerical / circumstellar matter
© The Authors 2024
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
The study of the continuum radiative transfer problem is crucial for the characterisation of circumstellar environments. Radiative processes play a major role in the determination of physical observables such as, for example, the temperature, abundances and velocity fields. The description of the radiation field is both a theoretical and numerical challenge, especially in the frequency-dependent and multi-dimensional case.
Several directions have been followed to tackle this problem. One approach involves approximate methods. It is usually done by assuming a particular form for the radiation field, most of the time based on symmetry arguments. The radiative transfer equation is then recast into a presumably simpler equation. This is the case, for example, for the Eddington approximation, in which the radiation is assumed to be a correction to an isotropic field, yielding an accurate description in the optically thick regime. For this approximation, the radiation is described by a simple linear diffusion equation. More sophisticated approximations were later developed, for example the flux-limited diffusion (Levermore & Pomraning 1981), which is asymptotically correct to both optically thin and thick regimes. In the latter case, the radiation is described by a non-linear diffusion equation.
While approximate methods are easier to handle numerically, they often fail to describe the radiation field accurately inside complex geometries (Kuiper & Klessen 2013). For these cases, it is necessary to numerically solve the radiative transfer equation directly (see Steinacker et al. 2013, for a thorough review of the different methods). A first approach is to approximate the transport operator with finite difference (e.g Steinacker et al. 2003), yielding a system of linear equations. This technique has the disadvantage of introducing spurious numerical oscillations and possible negative values for the intensity, due to the strong spatial and angular variations of radiation field. Other techniques, such as long and short characteristics (e.g Woitke et al. 2009), rely on the integral form of the radiative transfer equation. They are generally harder to implement than finite difference methods, are numerically demanding, and can also exhibit pathological behaviours such as negative values (Kunasz & Auer 1988). However, both approaches allow explicit error control. Finally, Monte-Carlo methods (e.g Wolf 2003) are amongst the most popular ones because they are easy and fast to implement, can handle complex geometries, and do not suffer from the same flaws as the other methods. However, they introduce noises that are inherent to their statistical nature. This class of methods additionally suffers from a reliable way of estimating the error, especially in optically thick regions where radiation is trapped inside the medium.
Aside from these techniques, finite element methods have already been used to solve the radiative transfer equation. A variant of it, the discontinuous Galerkin finite element method (Reed & Hill 1973, DGFEM hereafter), makes use of discontinuous elements and flux integrals along their boundaries, ensuring local flux conservation. However, as opposed to the classical finite element methods, the reconstructed solution is discontinuous across element edges. One of the main strengths of the method is its ability to produce a high-order numerical scheme, meaning a coarse computational grid can be used to achieve a small error. This feature is particularly interesting in the context of radiative transfer, where we are often limited by computational resources. Furthermore, the solution presents little to no oscillations, even in the cases where the solution displays a non-smooth behaviour (Nair et al. 2011), which often happens in the case of radiative transfer.
The DGFEM was successfully applied to the one-dimensional (1D) spherical transport problem, in the context of neutron transport (Machorro 2007), or more recently in grey stellar atmospheres (Kitzmann et al. 2016). Extensions to two-dimensional (2D) radiative transfer have been developed, in cylindrical coordinates, in the context of a non-local thermodynamic equilibrium (non-LTE) multilevel atom. Dykema et al. (1996) (and references therein) used a DGFEM with spatial linear elements for the evaluation of the Eddington tensor. The Eddington tensor was then subsequently used for closing the radiation field tensor moment system of equations. In another study, Cui & Li (2005) developed a DGFEM with a linear basis function in space and a step function in angle, for the computation of the radiative transfer in participating media.
In this study, we present the DGFEM applied to the frequency-dependent continuum radiative transfer problem, with isotropic scattering and coupled with the radiative equilibrium equation. We then used this method in the context of axis-symmetric dusty circumstellar environments. We show that this method can successfully determine the correct temperature profile, and allows for accurate images and spectral energy distributions (SEDs) to be computed by subsequent ray-tracing techniques. The paper is organised as follows: in Sect. 2 we describe the radiative transfer problem for axis-symmetric configurations. In Sect. 3 we present the DGFEM, and some of its numerical implementation features in Sect. 4. In Sect. 5, we compare this method against spherically and axis-symmetric benchmarks from the literature. Finally, in Sect. 6, we conclude and present some perspectives.
2 Description of the problem
We want to describe the radiation field inside an axis-symmetric circumstellar envelope. A central star of radius R* lies inside an inner cavity, free of matter, with a radius of Rin. The envelope spans from Rin up to the outer radius Rout. We assume that no matter is present after this radius. We consider the material to be exclusively made of dust. We choose to describe the problem with the spherical coordinate system, both for the spatial and angular variables (see Fig. 1). In addition to the axial symmetry around the polar axis, we also assume a planar symmetry with respect to the equatorial plane, at Θ = π/2. Given the symmetries, the domain of definition of radiation is D ⊂ ℝ4, with D ϶ x = [r, Θ, μ = cos θ, φ] ∈ [Rin, Rout] × (0, π/2] × [−1, 1] × [0, π]1. The radiation field is described by the time-independent radiation transfer equation (for a thorough derivation see e.g. Pomraning 1973, II-5), (1)
The subscript v denotes the frequency dependence, Iv is the specific intensity, and Ω.∇ is the transport operator, which corresponds to a spatial derivative in the direction Ω. We note that is the extinction coefficient, with and being the absorption and scattering coefficients, respectively. For dust species, we generally express these coefficients as , where are the absorption and scattering optical cross-sections and n is the number density. Dust grains are usually described as homogeneous spheres and their optical cross-sections are commonly computed with the help of Mie theory (Bohren & Huffman 1998). The emissivity ηv includes the thermal emission (Kirchoff-Planck law), proportional to the Planck function, Bv, at the local temperature T, and a scattering term. For this study, we assume the scattering to be isotropic, for the purpose of our numerical tests (although it can be generalised to anisotropic scattering with no difficulty). The emissivity is then , with Jv, the mean specific intensity, which is the zeroth-order angular moment of the specific intensity Iv, (2)
Additionally, the radiation field and the dust temperature are coupled via the equation of radiative equilibrium, (3)
The circumstellar matter is illuminated by a central star and it is customary to decompose the radiation field into two contributions (see e.g Steinacker et al. 2003), , where is the direct unprocessed stellar radiation field attenuated by the circumstellar extinction, and is the radiation emitted by the envelope (either via thermal emission or scattering). Eq. (1) can then be recast into the following system of equations, (4)
with . We note that in Steinacker et al. (2003), the thermal emission term was put on the right-hand side of the first Eq. (4). However, having it in the second equation is advantageous if , and are independent of temperature. In such cases, the equations decouple and can be solved and computed definitively. The first Eq. (4) can be integrated along a given ray, yielding the formal solution for the stellar contribution, (5)
The incident stellar radiation field is . Again, for the purpose of the numerical tests, we assume the star to radiate as a black body at the temperature T*. Furthermore, μ* = (1 − (R*/r)2)1/2 is the cosine of the angle subtended by the star at radius r, and s(Rin) is the distance between a given point r and Rin, in the direction Ω. The argument in the exponential is the opposite of the optical depth integrated along the ray. For dusty media, we usually have r ≫ R*, thus the star can be treated as a point source. In the point source approximation, μ* ≈ 1 − (R*/r)2 /2, and the stellar mean intensity can be expressed analytically as, (6)
In general, the integral in Eq. (6) can be carried out numerically, providing the stellar source term for Eq. (4).
To complete the description of the problem, we specify the boundary conditions for . We do so by prescribing the incident intensity upon the surface of the domain D located at rs, (7)
where ŝ is the unit vector normal to the surface of D, pointing outside of the domain. At the inner radius, , and the incident radiation, in a given direction , comes directly from the opposite point of the cavity, (8)
with (Rin, Θ′, μ′, φ′) being the coordinates of the opposite point. Their derivation is given in Appendix A. On the outer edge, , and we assume that there is no incident radiation upon the surface, (9)
At the equator, , and the planar symmetry requires the radiation field to be, as shown in Fig. 2, (10)
Fig. 1 Illustration showing the coordinate system. We note that is the standard spherical basis. Given the symmetry around the polar axis, the radiation field, at a given position r, in a given direction Ω is a function of two spatial (r, Θ) and two angular (μ = cos θ, φ) coordinates. |
Fig. 2 Representation of the direction vector Ω at the equator. The dotted arrows represent the symmetric of Ω with respect to the equatorial plane . The left and right 2D views allow the corresponding direction angles to be seen. |
3 The radiative transfer equation with the discontinuous Galerkin finite element method
We present the DGFEM applied to the radiative transfer Eq. (1). We used some elements of notation from Kitzmann et al. (2016) and Hesthaven & Warburton (2007). For simplicity, we omitted the frequency subscript. The conservative form of Eq. (1) is, (11)
with . We introduced the variable , which is the quantity that is being conserved in Eq. (11). Using this quantity is important because it improves the stability of the numerical scheme, especially near the polar axis (Θ = 0) where Eq. (1) is not defined when using the spherical coordinate system. The transport operator, ∇x., corresponds to the Cartesian divergence operator with respect to x = (r, Θ, μ, φ), applied to the flux vector F of , (12)
We decomposed the domain D into N = Nr × NΘ × Nμ × Nφ non-overlapping rectangular elements Di,j,k,l, with Nr, NΘ, Nμ and Nφ being the number of elements along the r, Θ, μ and φ coordinates, respectively. Each element is denoted with the help of four indexes i, j k, and l, ranging from 0 to Nr − 1, NΘ − 1, Nμ − 1 and, Nφ − 1. Inside each element Di,j,k,l, we used the nodal representation and we approximated the local solution by a four-dimensional (4D) polynomial expansion, (13)
with na, nb, nc, and nd, being the number of nodes inside each element Di,j,k,l, along the r, Θ, μ, and φ coordinate, respectively. We note that ha,b,c,d is the 4D Lagrange polynomial, defined as, (14)
By definition, the coefficients correspond to the value of at the nodes of coordinates xa,b,c,d = (ra, Θb, μc, φd). An example of an element Di,j,k,l with the associated nodes is shown in Fig. B.1. The global numerical approximation of the solution Ih across the domain D is formed by the sum of the N piece-wise continuous solutions inside each element, (15)
Now we shall introduce ℛh, the residual of Eq. (11), (16)
with . We form the classical Galerkin formulation (see e.g. Eq. (2.3) of Hesthaven & Warburton 2007) by requiring the residual to be orthogonal, inside each element, to the same set of functions used for the solution representation, (17)
The divergence term that appears in Eq. (17) can be recast with the help of the divergence theorem (Green-Ostrogradsky), yielding the following so-called weak formulation of the radiative transfer equation Eq. (11): (18)
The first term, in the left-hand side of Eq. (18), is a surface integral, along the boundaries of Di,j,k,l. Furthermore, ŝ is the outward normal vector to the surface element, and F* is an estimate of the flux at the cell interface, referred to as the numerical flux. It arises because the solution is not uniquely defined at the edge of the element, due to the discontinuous nature of the solution. We would like to emphasise that Di,j,k,l is a 4D rectangular element, and hence that each element has 2 × 4 boundary surfaces. The second and third terms in Eq. (18) are volume integrals and are purely local terms. Consequently, they only depend on the solution inside the element considered. For F*, several choices are possible, depending on the nature of the problem (Cockburn 2003). In our case, and as is commonly employed for transport problems, we used the Lax-Friedrichs numerical flux, defined as (e.g. on the radial right element edge where ), (19)
where and denote the left and right values of at the element edge, respectively, and ar and Fr are as defined in Eq. (12). The same form of expression holds for the other surfaces of the element.
Eq. (18) can be assembled into a system of N′ = N × n equations with n = na × nb × nc × nd, relating the coefficients. The integrals are numerically estimated with the help of a quadrature formula. The choice of the quadrature, with the associated roots (nodes), is not unique and is usually problem-dependent (for an extensive review see e.g Kopriva & Gassner 2010). In general, we can put the system of Eq. (18) in the form of the following linear system: (20)
with being the vector of size N′, which contains the solution points for the full domain D, , a sparse matrix with a size of N′ × N′ coupling the elements of , and ℬ a vector of size N′ containing the emissivity term .
4 Solution strategy and numerical considerations
The solution of the problem presented in Sec. 2 involves determining the radiation field, given by Eq. (4), coupled with the equation of radiative equilibrium, Eq. (3). For the test cases we consider in Sec. 5, the star was treated in the point source approximation, and hence we directly used Eq. (6) to compute . The radiation field of the envelope , described by the second Eq. (4), was solved with the DGFEM presented in Sec. 3. For simplicity, we omitted the envelope superscript but implicitly refer to this contribution below.
Because the radiative transfer problem is intrinsically highly dimensional, the size of the matrix is often huge and solving Eq. (20) becomes numerically tedious. However, is sparse, because each solution point, , inside a given element, Di,j,kl, only depends on the other points inside the same element and the neighbouring ones, inside Di±1,j±1,k−1,l+1. This property allows us to rewrite Eq. (20) as, (21)
We note that i,j,k,l are the diagonal blocks of , with a size of n × n while i±1,j±1,k−1,l+1 (with a size of n × n) are the only nonzero, non-diagonal blocks. Furthermore, the elements Di,j,k+1,1 and Di,j,k,l−1 do not contribute because of the expression of the Lax-Friedrichs numerical flux, Eq. (19), with aμ ≥ 0 and aφ ≤ 0, ∀x ∈ D. We note that and ℬi,j,k,l are the sub-vector of and ℬ, respectively, with a length of n, containing the local points in Di,j,k,l.
To put in a matrix form, we used a global index α = i NΘ Nμ + j Nμ Nφ, + k Nφ + l. An example of the structure of the matrix with the associated blocks is displayed in Fig. 3. The formulation Eq. (21) avoids the storage and computation of the entire matrix , which reduces the computational effort.
In general, the simplest approach to obtain a solution from the problem is to solve Eq. (21) and to update the temperature with Eq. (3). Iterating between these two steps until convergence yields the solution to the problem. This procedure, commonly referred to as the Λ-iteration in the literature, becomes very slow and does not converge for large optical depths (Mihalas & Weibel-Mihalas 1999, VI-83). Our solution strategy is directly inspired from Perdigon et al. (2021). The key point of the method is to solve simultaneously instead of repetitively, Eq. (21) and Eq. (3). This strategy can be assimilated to an acceleration procedure to the usual Λ-iteration and yields satisfying results up to moderately thick envelopes. We however note that, as for the usual Λ-iteration, it converges very slowly for optically thick envelopes.
We proceeded with the following solution strategy: if we denote the iteration of the method by the superscript index n, we first (i) computed the stellar mean radiation field with Eq. (6) and set , as an initial condition. This allowed us to compute the initial temperature profile T0, with the help of Eq. (3), where . Then (ii), for each frequency, we computed with the help of Eq. (21), which we rewrote, performing a block Gauss-Seidel sweep (see e.g Karniadakis & Kirby II 2003, VII-2), (22)
We give in appendix B the expressions of i,j,k,l and of the right-hand side of Eq. (22). We note that Eq. (22) represents N linear systems to solve. We solved the linear system Eq. (22) by direct inversion, using the Gauss elimination algorithm (see e.g Karniadakis & Kirby II 2003, IX-1). (iii) We computed and consequently updated the temperature via Eq. (3). The new temperature allowed us to update the right-hand side of Eq. (22) and to repeat steps (ii) and (iii), until convergence. We would like to point out that this scheme is not strictly identical to the usual Λ-iteration as we updated the temperature together with Iv after each iteration index n. Finally, concerning the implementation of the boundary conditions, we directly followed the prescription from Kitzmann et al. (2016, Appendix A.6).
Summary of the main results.
Fig. 3 Example of the sparse structure of with Nr = NΘ = Nμ = Nφ = 2. The squares represent the blocks of size n × n containing nonzero values. The matrix has bands corresponding to the coupling blocks , and . |
5 Numerical tests
The lack of analytic solutions for the radiative transfer problem, especially in the multi-dimensional and frequency-dependent case, limited our options for testing the validity of the method. In general, one has to compare numerical solutions with other ones from the literature. For this purpose, benchmark problems have been proposed. The first test case we consider, in Sect. 5.1, is the frequency-dependent radiative transfer problem inside a spherically symmetric envelope, from Ivezic et al. (1997). The second test case, in Sec. 5.2, is about the frequency-dependent radiative transfer inside an axis-symmetric envelope (disc), from Pascucci et al. (2004). For the latter test, we note that five different codes were used (including Monte Carlo and grid-based methods) to produce the benchmark and check the consistency of the results. We note that both of these tests are compatible with the boundary conditions presented in Sec. 2. A summary of the main results is presented in Table 1.
5.1 1D spherically symmetric envelope
A point source, surrounded by a spherically symmetric envelope of dust, at radiative equilibrium, radiates as a black body at the temperature T* = 2500 K. This envelope extends from the inner radius Rin to the outer radius Rout = 1000 Rin. The inner radius was set so that its temperature would always be Tin = T(Rin) = 800 K. The number density profile, n(r), is assumed to be a power law, of the form n(r) = n0 (Rin/r)2. The number density at the base of the disc, n0, is determined by setting the value of , which corresponds to the radial optical depth integrated though the envelope, at v0 = c/λ0, with λ0 = 1 μm, (23)
We note that is the extinction cross-section coefficient at v0. In this test we consider and 102 to correspond to a moderately thin and thick envelopes, respectively. The absorption and scattering cross-sections, and , feature a bi-linear behaviour in log-log scaling, with a constant profile for v ≥ v0 and a power-law dependence ∝ (v/v0)α, for v ≤ v0, with α = 1 and 4 for absorption and scattering, respectively. This dependence aims to mimic the behaviour of spherical astronomical silicate grains. We set the value of the thermal coupling parameter to be 1/2. The benchmarks from Ivezic et al. (1997) were reproduced with version 2 of DUSTY2. Although the envelope is spherically symmetric, we used a grid of N = Nr × NΘ × Nμ × Nφ = 164 elements in our DGFEM code, with each elements containing na ×nb × nc × nd = 3 × 2 × 3 × 3 nodes, as pictured in Fig. B.1. For the radial coordinate, the cell edges were logarithmically spaced, to account for the important dynamic of the solution with respect to the radius. The frequency grid consists of 60 logaritmically spaced points, ranging from λ = 10−2 μm to λ = 3.6 × 104 μm.
The temperature profiles, T, and the SEDs, λFλ/F with , of the envelope are shown in Fig. 4. For the DGFEM code, the temperature that is shown corresponds to the mean radial profile across all angular points Θ. Additionally, the normalised SEDs were computed with the help of a ray-tracing module we present in Appendix C. We subsequently explain the reasons for such a choice in Sect. 5.3. The spatial and frequency grids differ between both codes and we performed a linear interpolation (in log-log scaling) of the DUSTY profiles at our grid points, in order to do the comparison. We observed a good agreement between the two codes. On average, the absolute relative differences stay below 1% for the temperatures and 2% for the SEDs.
Fig. 4 Temperature profiles (left panels) and normalised SEDs (right panels) for the spherically symmetric envelope, with (blue curve) and (orange curve). The cross marks represent the solution from this study and the dashed curves are from DUSTY (Ivezic et al. 1997). The lower panels show the relative differences between the two codes. |
5.2 2D axis-symmetric envelope
A point source surrounded by an axis-symmetric disc of dust, at radiative equilibrium, radiates as a black body at the temperature T* = 5800 K. This disc extends from the inner radius Rin = 1 AU to the outer radius Rout = 1000 AU. The density profile n(r, Θ) is assumed to be (24)
with , and zd = Rout/8. This density law is characteristic of a Keplerian disc hydrostatically supported in the vertical direction, assuming a constant temperature along this direction (Chiang & Goldreich 1997). Again, n0 was determined to set , the radial optical depth, integrated through the disc mid-plane (Θ = π/2), at frequency v0 = c/λ0 with λ0 = 0.55 μm, (25)
The opacities were taken from Draine & Lee (1984). They are associated with a unique dust species composed of spherical astronomical silicate grains, with a radius of 0.12 μm and a density of 3.6 g cm−3. A table of pre-computed values for and is available in Pascucci et al. (2004). We considered the optically thin and thick cases and 102, respectively. The benchmarks were produced with RADMC-3D3 (Dullemond et al. 2012), a Monte-Carlo radiative transfer code and whose previous version is part of the original benchmarks from Pascucci et al. (2004). For our code, we used the same spatial and frequency grids as in the previous test. For RADMC-3D, we used 128 points both in radial and angular direction, with a logarithmic sampling in the radial direction.
The temperature of the disc is displayed in Fig. 5. The DGFEM code successfully reproduces it. The regions of the disc with steep gradients are always well reproduced, even with a fairly reasonable number of nodes. This result is a direct consequence of having a high-order numerical scheme. The radiation field can exhibit discontinuities, because of boundary conditions or very strong density gradients. Our numerical tests revealed little to no oscillations and very few negative values for the specific intensity, which did not pollute the computation of the mean radiation field. On average, the absolute relative differences stay below one percent for both test cases. The temperature in the disc mid-plane is very well reproduced, highlighting that the method is able to correctly represent the shadowed regions of the disc (the cold outer mid-plane regions shadowed by the inner dense regions).
In Fig. 6, we show the corresponding emerging SEDs, for two inclinations with respect to the polar axis i = 12.5 and 77.5 deg. Again our SEDs were computed from the ray-tracing of the emissivity ηv. They consist of a stellar and an envelope contribution (see Eq. (C.1)). The stellar component that peaks at around λ ≈ 0.6 μm was computed from of Eq. (C.2) and the discrepancies are always < 1% where this part dominates. Concerning the contribution from the envelope, both the scattering (< 1 μm) and emission (>10 μm) parts agree well. On average, the absolute relative differences stay below 5%, which are in the typical discrepancy levels of the original benchmark. We note a peak in the discrepancies for the optically thick case, between 8.7 and 10.2 μm (according to the resolution of our frequency grid). The same behaviour was previously observed in Pascucci et al. (2004), between a grid-based and a Monte-Carlo code. The authors suggested that, at these wavelengths, the flux mainly comes from the inner disc regions (between 1 and 2 AU) and the numerical simulations are particularly sensitive to the resolution of the inner parts. However, we tried to increase the resolution in these regions, which did not result in any improvement.
In Fig. 7, we show a set of images from the DGFEM code of the 10 AU disc inner regions, at λ = 2.3, 4.5, and 12.1 μm and for several inclinations i = 12.5, 77.5, and 90 deg. These wavelengths are characteristic of the operating spectral bands of the interferometric instruments, such as GRAVITY (Gravity Collaboration 2017) and MATISSE (Lopez et al. 2022). On average, the agreement between the images from the DGFEM code and the images from RADMC-3D is around 10%, for all frequencies and inclinations. We show in Fig. 8 the comparison of an image slice at x = 0 AU, for two inclinations i = 12.5 and 77.5 deg and for λ = 2.3, 4.5, and 12.1 μm. In general, the disc emitting inner regions (peaks in Fig. 7) are reproduced very well . The biggest discrepancies occur in the wings of the peaks, where the gradient of intensity is the steepest (in logarithmic scaling).
Fig. 5 Temperature profiles for the axis-symmetric envelope with (blue curve) and (orange curve), in the in the disc midplane (left panel) and at r = 2 AU (right panel). The cross marks represent the solution from this study and the dashed curves were computed with RADMC-3D. The lower panels show the relative differences between the two codes. |
Fig. 6 SED profiles for the axis-symmetric envelope with (blue curve) and (orange curve). The left and right panels correspond to i = 12.5 and 77.5 deg, respectively. The cross marks represent the solution from this study and the dashed curves were computed with RADMC-3D. The lower panels show the relative differences between the two codes. |
Fig. 7 Images at λ = 2.3, 4.5, and 12.1 μm (left, middle, and right panels, respectively) of the 10 AU inner regions of the axis-symmetric envelope (Sect. 5.2) and computed from the DGFEM solution with the ray-tracing module (Appendix C). The top, middle and bottom panels correspond to the inclinations i = 12.5, 77.5, and 90 deg, respectively. The colour code shows the specific intensity value (in W m−2 Hz−1 sr−1) of the envelope , inside each pixel. The solid white lines show the iso-contours of . |
5.3 Effect of the grid resolution on temperatures and SEDs
As we previously mentioned in Sects. 5.1 and 5.2, we used a ray-tracing procedure to accurately compute SEDs. However, in practice, we already know the intensity Iv(r = Rout, Θ, μ, φ) from the DGFEM solution at the outer edge of the envelope. This intensity can theoretically be used to compute the SEDs directly with no further post-processing treatment required. It is, however, not advised to proceed in this way, at least in multi-dimensional cases where the radiation field has a complex geometry, because of the poor accuracy of the local solution on the outer edge due to the poor resolution of the grid (Rout ΔΘ is large). This effect becomes even more important if we are interested in studying the inner parts of the disc and hence probing a very small part of the outer shell of the envelope. Consequently a post-processing ray-tracing step is usually needed in order to capture all the features of the disc accurately. The post-processing ray-tracing step requires a good estimate for the emissivity (or source function, see Eq. (11) in the whole envelope to compute accurate SEDs (see Appendix C). The emissivity is implicitly a function of the temperature of the medium (through the black-body emission Bv(T) and the radiative equilibrium Eq. (3)), and hence an accurate temperature is needed to compute SEDs with precision.
In order to study the effect of the resolution of the computational grid on the temperature agreement, we computed the L2 norm of the relative residual of the temperature, between our code TDG and the RADMC-3D code Tref. It is defined as, (26)
where the integral runs across the volume of the envelope. In Fig. 9, we display the L2 norm versus the total number of elements of our grid (we would like to emphasise that the total number of elements is N = Nr × NΘ × Nμ × Nφ and that each element has 3 × 2 × 3 × 3 nodes). We varied N in such a way as to keep the number of elements along each dimension identical (Nr = NΘ = Nμ = Nφ).
We can see that the L2 norm saturates when N ≈ 104, which corresponds to a number of elements along each dimension of ten. After that point, no substantial improvement in the temperature agreement is obtained. The computational grid at this point is still quite coarse, compared to the commonly used methods in the literature such as finite differences (Steinacker et al. 2003), short-characteristics, (Dullemond & Turolla 2000) or Monte-Carlo (Wolf 2003) methods. This is one of the advantages of the DGFEM that can achieve a given error threshold with less resolution, thanks to the form of the intensity inside each element (see Eq. (13)).
In Fig. 10, we compare the SEDs obtained either from the post-processing ray-tracing or directly by using Iv(r = Rout, Θ, μ, φ). We compare these SEDs with the benchmark curves we already presented in Sects. 5.1 and 5.2 (τ = 100). While for the spherically symmetric case, they agree quite well, Iv(r = Rout, Θ, μ, φ) fails to give an accurate SED for the 2D disc configuration. This is expected since for this case, the radiation field has a dependence on the Θ and φ variables. We note that the increase in the grid size should make the SED from Iv(r = Rout, Θ, μ, φ) converge to the ray-tracing one, but we could not investigate this hypothesis due to the limits of our current computational resources. In the absence of a finer grid, the post-treatment ray-tracing is thus needed if one wants to compute the emerging fluxes from the disc with precision.
Fig. 8 Image slices at x = 0 AU, for two inclinations i = 12.5 and 77.5 deg (left and right panels, respectively). The blue, orange, and green curves are the RADMC-3D intensities at λ = 2.3, 4.5, and 12.1 μm, respectively. The cross marks represent the solution from this study. The lower panels show the relative differences between the two codes. |
Fig. 9 L2 norm of the relative differences in temperature between this study and RADMC-3D versus the number of elements in the computational grid (N = Nr × NΘ × Nμ × Nφ. The L2 norm was computed for the benchmark problem presented in Sec. 5.2. |
6 Conclusions
For this study, we applied the DGFEM to the frequency-dependent 2D radiative transfer equation, coupled with the radiative equilibrium equation, inside axis-symmetric circumstellar envelopes. We have shown that it can accurately compute the temperature field and allow for the correct determination of images and SED profiles, via ray-tracing techniques. A desirable feature of the method is the ability to control the order of the numerical scheme via the number of nodes inside each element, meaning that a high-order numerical scheme can simply be achieved by increasing the number of nodes inside each element.
The DGFEM formulation, Eq. (18), is particularly adapted to parallelisation and adaptive mesh refinement (AMR) grids (e.g. Frisken & Perry 2002). In this regard, direct access to the residual, Eq. (16), might provide a robust estimate for the error that could be used as a criterion for the grid refinement (see also Richling et al. 2001). We also note that a further 3D generalisation is also possible and straightforward. Together with Monte-Carlo and ray-tracing techniques, the DGFEM provides an additional method for solving the radiative transfer problem, and it could be used in cases where the other methods are expected to be less efficient.
Fig. 10 Emerging SEDs of RADMC-3D (solid black lines) and of the DGFEM code, computed from the post-processing ray-tracing procedure (orange dots), and from the DGFEM solution at the outer edge (blue crosses). The left panel corresponds to the optically thick spherically symmetric case (τ = 100, see Sect. 5.1) while the right panel corresponds to the optically thick disc benchmark from Sect. 5.2. The lower panels show the relative difference between the curves of this study and the benchmarks. |
Acknowledgements
This study has been supported by the Lagrange laboratory of Astrophysics and funded under a 3-year PhD grant from Ecole Doctorale Sciences Fondamentales et Appliquées (EDSFA) of the Université Côte-d’Azur (UCA). The 1D benchmark profiles were computed with the radiative transfer code DUSTY, available at http://faculty.washington.edu/ivezic/dusty_web/. The 2D benchmark were computed with RADMC-3D: https://www.ita.uni-heidelberg.de/~dullemond/software/radmc-3d/ Part of the computations were carried out with the help of OPAL-Meso computing facilities. The authors are grateful to the OPAL infrastructure from Observatoire de la Côte d’Azur (CRIMSON) for providing resources and support. We are grateful to the referee and the editor whose insights helped improve this work. We also want to express our gratitude to the language editor who helped improve the clarity of this paper.
Appendix A Boundary conditions on a spherical enclosed cavity
The boundary condition in Eq. (7) requires to know, for a given point A (Rin, Θ) on the inner cavity of radius Rin and in a direction , the coordinates of the opposite point B (Rin, Θ′, μ′ = cos θ′, φ′). In the following, the coordinates with a prime superscript are associated with the opposite point B. The problem is illustrated in Fig. A.1, where we show the plane containing the point A, B and the star. The problem is axis-symmetric around the polar axis and the radiation field does not depend on the azimutal angle Φ (see Fig. 1). Consequently, we can arbitrarily choose the Φ coordinate of the point A and set it to be Φ = π/2. This convention simplifies the computations.
Fig. A.1 Geometry of a ray in the inner cavity. Point B is the opposite point of A, in the direction Ω. |
First, we see that θ′ = π − θ because the OAB triangle is isosceles, we then have, (A.1)
The position rB of the opposite point B, is linked to the position rA of the point A through the relation (A.2)
with sAB the distance between the point A and B. This length is equal to sAB = 2 μ Rin, because it corresponds to the base of the isosceles triangle OAB, with OA = OB = Rin. The direction vector, Ω, is, (A.3)
Eq. (A.2) is then (A.4)
From the z-component of Eq. (A.4) we get, (A.5)
and consequently Θ′ = arccos | cos Θ′|. The absolute value occurs if we consider the planar symmetry with respect to the equatorial plane.
The angle φ′ at the point B verifies, (A.6)
and Ω given by Eq. (A.3). Then, using Eq. (A.4) in Eq. (A.6) yields, after tedious calculations, (A.8)
Consequently φ′ = arctan |u|/υ with u and υ being the numerator and denominator in Eq. (A.8), respectively. This time, the absolute value is present because the symmetry around the polar axis causes the φ′ angle to take its values in [0, π]. We note that the previous formula is not valid for the couples of coordinates (μ, φ) = (sin (Θ/2), 0) and (μ, φ) = (cos (Θ/2), π), because it corresponds to have the opposite point B on the poles, where and are not defined. These particular cases are in practice not a problem since the radiation field is independent of φ′ at these points.
Appendix B DGFEM calculations
The computation of the terms in Eq. (18) requires the choice of a particular quadrature in the cell Di,j,k,l, for each coordinate (r, Θ, μ, φ). For the r, μ, φ coordinates, we make use of the Gauss-Lobatto quadrature, which has the advantage of having roots at the end points of the cell. This avoids the use of an interpolation formula when computing the numerical flux at the element edges. The Gauss-Lobatto quadrature can exactly integrate polynomials up to degree 2 n − 3, with n the number of nodes. For the coordinate Θ, we cannot use the Gauss-Lobatto quadrature, at least in the cells that are touching Θ = 0, because the radiative transfer equation in the spherical coordinates system is not defined at the pole. For this reason and to keep an homogeneous method along Θ, we use a Gauss-Legendre quadrature for this coordinate. The Gauss-Legendre quadrature is more precise and can exactly integrate polynomials up to degree 2n − 1, but at the expense of an interpolation method required to compute the flux at the cell edges. An example of an element Di,j,k,l with the chosen nodes is shown in Fig. B.1.
Fig. B.1 Example of an element Di,j,k,l (dashed lines) using na = nc = nd = 3 and nb = 2 points. The left and right panels correspond to 2D slices in the (r, Θ) and (μ, φ) planes, respectively (we display here the local coordinate system given by Eq. B.2). The black dots correspond to the nodes where the solution is computed while the crosses correspond to the interpolated value of I used to compute the numerical flux at the interface, along the Θ coordinate. |
In the following, all the superscript indexes refer to the element identification while the subscripts denote each node in the considered element. We start with the volume term in Eq. (18), (B.1)
We note that Δxi,j,k,l = Δri, Δ Θj Δμk Δφl is the 4D volume of the element Di,j,k,l For integration, we rather use the local coordinates , defined as (e.g for the r coordinate), (B.2)
with Δri, the element width along the coordinate r and ri+1/2, the radial coordinate of the centre of the element. The same expression holds for the other coordinates. The quantities are the weights associated with the different quadrature in each direction. Finally, is the 4D Lagrange polynomials, defined in Eq. (14), evaluated at the node . By definition of the Lagrange polynomials, we have, (B.3)
with δa′,a the usual delta Kronecker. The other volume term in Eq. (18) is expressed as, (B.4)
In Eq. (B.4), we used the definition of the flux Eq. (12). The quantity is the partial derivative of the Lagrange polynomial, with respect to the coordinate r, evaluated at the node (with similar definitions for the other coordinates).
The last term to evaluate is the surface integral (first term in Eq. 18). The 4D element is delimited by 2 × 4 = 8 surfaces. We give here the derivation for the surface integrals normal to the coordinate r and the other terms will follow by substitution of the indices. We have for the radial right and left surfaces , respectively, (B.5)
At the cell boundaries, we use the upwind numerical flux, for example, at the right edge, (B.6)
The use of the Gauss-Lobatto quadrature further simplify the computation because the quadrature nodes include the endpoints . The numerical flux can then be expressed as (B.7)
The same form holds for the surface integrals normal to the μ and φ coordinates. We however note that we have aμ ≥ 0, aφ ≤ 0 ∀ x ∈ D, so there are no terms proportional to and . For the flux computation normal to the Θ coordinate, we do not directly have the value of the solution at the element interface (see Fig. B.1), we need to interpolate the solution with the help of Eq. (13), (B.8)
All the terms in this section can be put in the form of the system of equations, given by Eq. (22), (B.9)
In Eq. (B.9), we replaced the RHS by bi,j,k,l because we do not need to formally write the non-diagonal matrices i±1,j±1,k−1,l+1. To assemble i,j,k,l, we make use of the global index α = a nbncnd + b ncnd + c nd + d. The elements of i,j,k,l and bi,j,k,l are then, (B.10) (B.11)
Appendix C Ray-tracing module
The SEDs and intensity maps from the DGFEM code, shown in Sect 5, were computed with the help of the ray-tracing routine we present here. This procedure is quite generally used in the literature, for all type of codes, and is not a limitation of the DGFEM method itself (see e.g Pinte et al. 2009, Sect 2.2.3)
For the SED, we need to estimate the total flux that an observer receive from the object, situated at a distance d ≫ Rout and doing an angle i with the polar axis (see Fig. C.1). Because we decoupled the stellar from the envelope radiation (see Eq. 4), the total flux is made of the stellar and envelope total flux, (C.1)
If we assume the star to be an unresolved black-body point source, the stellar flux at distance d is the flux of the star at the stellar surface attenuated by the circumstellar matter present in the direction of the line of sight (black dotted line in Fig. C.1), and with the dilution factor (R*/d)2, (C.2)
Fig. C.1 Example of a ray (red line) normal to the image plane, crossing the spherical grid. In this example, we display a square image of size L × L and a circular image of radius R. The red crosses represent the intersections between the ray and the grid. The values of and Sv at these intersections are linearly interpolated from the grid adjacent values (red dots). |
To compute the envelope flux and intensity maps, we define an image plane , at a distance d ≫ Rout from the star and tilted with an angle i with respect to the polar axis. The x and y axes are oriented with the help of the spherical coordinates system . In this plane, we can construct images of any geometry but we only consider the special cases of a square image of width L and a circular image of radius R. Since d ≫ Rout, the flux inside the image can formally be written, (C.3)
In practice, images are made of a collection of pixels in which we evaluate the emerging specific intensity at the pixel centre. A square (circular) image, is divided into N × N (Nr × Nω) pixels and the flux in the image can then be rewritten, (C.4)
with Δxi Δyi (ri Δwj Δri;) the pixel size of the square (circular) image. We note that the circular image is particularly well-suited for the computation of since we can easily increase the number of pixels in the centre of the image in order to resolve the disc inner parts.
The emerging specific intensity crossing each pixel centre for a circular image), along the ray normal to the image plane (red line in Fig. C.1) is computed by integration of the emissivity along the ray, (C.5)
We define s as to be the distance from the pixel centre r0 to a given point along the ray. The quantities ηv and are the emissivity and the extinction coefficient, respectively, as defined in Eq. (1). The points smin and smax correspond to the two intersections of the ray with the sphere of radius Rout.
In practice, ηv and are defined on a discrete grid and the previous integral can be rewritten as, (C.6)
The {si} are the n intersections between the ray and the grid (red crosses in Fig. C.1), with s0 = smin and sn−1 = smax. The coordinates of all intersections can be computed, in the Cartesian coordinate system. We can express , the contribution to the total intensity of each portion between two consecutive intersections as, (C.7)
Following Olson et al. (1986), we assume that ηv and are linear functions between two consecutive intersections. Each contribution is given by, (C.8)
and . The values of Sv and at the intersections si and si+1 are estimated by linear interpolation from the grid-adjacent values (red dots in Fig. C.1). The optical depth, , can be computed recursively, (C.9)
with and defined in Eq. (C.8).
References
- Bohren, C. F., & Huffman, D. R. 1998, Absorption and Scattering by a Sphere (John Wiley & Sons, Ltd), 82 [Google Scholar]
- Chiang, E. I., & Goldreich, P. 1997, ApJ, 490, 368 [Google Scholar]
- Cockburn, B. 2003, Z. Angew. Math. Mech., 83, 731 [NASA ADS] [CrossRef] [Google Scholar]
- Cui, X., & Li, B. Q. 2005, JQSRT, 96, 383 [NASA ADS] [CrossRef] [Google Scholar]
- Draine, B. T., & Lee, H. M. 1984, ApJ, 285, 89 [NASA ADS] [CrossRef] [Google Scholar]
- Dullemond, C. P., & Turolla, R. 2000, A&A, 360, 1187 [NASA ADS] [Google Scholar]
- Dullemond, C. P., Juhasz, A., Pohl, A., et al. 2012, RADMC-3D: a multipurpose radiative transfer tool, Astrophysics Source Code Library [record ascl:1202.015] [Google Scholar]
- Dykema, P. G., Klein, R. I., & Castor, J. I. 1996, ApJ, 457, 892 [NASA ADS] [CrossRef] [Google Scholar]
- Frisken, S., & Perry, R. 2002, J. Graph. Tools, 7, 1 [CrossRef] [Google Scholar]
- Gravity Collaboration (Abuter, R., et al.) 2017, A&A, 602, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Hesthaven, J. S., & Warburton, T. 2007, Nodal Discontinuous Galerkin Methods: Algorithms, Analysis, and Applications (Springer Science & Business Media), 19 [Google Scholar]
- Ivezic, Z., Groenewegen, M. A. T., Men’shchikov, A., & Szczerba, R. 1997, MNRAS, 291, 121 [NASA ADS] [CrossRef] [Google Scholar]
- Karniadakis, G. E., & Kirby II, R. M. 2003, Parallel Scientific Computing in C and MPI: A Seamless Approach to Parallel Algorithms and their Implementation (Cambridge University Press) [Google Scholar]
- Kitzmann, D., Bolte, J., & Patzer, A. B. C. 2016, A&A, 595, A90 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Kopriva, D. A., & Gassner, G. 2010, J. Sci. Comput., 44, 136 [CrossRef] [Google Scholar]
- Kuiper, R., & Klessen, R. S. 2013, A&A, 555, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Kunasz, P., & Auer, L. H. 1988, J. Quant. Spec. Radiat. Transf., 39, 67 [NASA ADS] [CrossRef] [Google Scholar]
- Levermore, C. D., & Pomraning, G. C. 1981, ApJ, 248, 321 [NASA ADS] [CrossRef] [Google Scholar]
- Lopez, B., Lagarde, S., Petrov, R. G., et al. 2022, A&A, 659, A192 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Machorro, E. 2007, J. Computat. Phys., 223, 67 [CrossRef] [Google Scholar]
- Mihalas, D., & Weibel-Mihalas, B. 1999, Foundations of Radiation Hydrodynamics, Dover Books on Physics (Dover Publications) [Google Scholar]
- Nair, R. D., Levy, M. N., & Lauritzen, P. H. 2011, Emerging Numerical Methods for Atmospheric Modeling (Berlin, Heidelberg: Springer Berlin Heidelberg), 251 [Google Scholar]
- Olson, G. L., Auer, L. H., & Buchler, J. R. 1986, J. Quant. Spec. Radiat. Transf., 35, 431 [NASA ADS] [CrossRef] [Google Scholar]
- Pascucci, I., Wolf, S., Steinacker, J., et al. 2004, A&A, 417, 793 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Perdigon, J., Niccolini, G., & Faurobert, M. 2021, A&A, 653, A139 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pinte, C., Harries, T. J., Min, M., et al. 2009, A&A, 498, 967 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Pomraning, G. C. 1973, The Equations of Radiation Hydrodynamics (Oxford: Pergamon Press) [Google Scholar]
- Reed, W. H., & Hill, T. R. 1973, Triangular mesh methods for the neutron transport equation, Tech. rep. (North Mexico, USA: Los Alamos Scientific Laboratory) [Google Scholar]
- Richling, S., Meinköhn, E., Kryzhevoi, N., & Kanschat, G. 2001, A&A, 380, 776 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Steinacker, J., Henning, T., Bacmann, A., & Semenov, D. 2003, A&A, 401, 405 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Steinacker, J., Baes, M., & Gordon, K. D. 2013, ARA&A, 51, 63 [CrossRef] [Google Scholar]
- Woitke, P., Kamp, I., & Thi, W. F. 2009, A&A, 501, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
- Wolf, S. 2003, Comput. Phys. Commun., 150, 99 [CrossRef] [Google Scholar]
Available at http://faculty.washington.edu/ivezic/dusty_web/
All Tables
All Figures
Fig. 1 Illustration showing the coordinate system. We note that is the standard spherical basis. Given the symmetry around the polar axis, the radiation field, at a given position r, in a given direction Ω is a function of two spatial (r, Θ) and two angular (μ = cos θ, φ) coordinates. |
|
In the text |
Fig. 2 Representation of the direction vector Ω at the equator. The dotted arrows represent the symmetric of Ω with respect to the equatorial plane . The left and right 2D views allow the corresponding direction angles to be seen. |
|
In the text |
Fig. 3 Example of the sparse structure of with Nr = NΘ = Nμ = Nφ = 2. The squares represent the blocks of size n × n containing nonzero values. The matrix has bands corresponding to the coupling blocks , and . |
|
In the text |
Fig. 4 Temperature profiles (left panels) and normalised SEDs (right panels) for the spherically symmetric envelope, with (blue curve) and (orange curve). The cross marks represent the solution from this study and the dashed curves are from DUSTY (Ivezic et al. 1997). The lower panels show the relative differences between the two codes. |
|
In the text |
Fig. 5 Temperature profiles for the axis-symmetric envelope with (blue curve) and (orange curve), in the in the disc midplane (left panel) and at r = 2 AU (right panel). The cross marks represent the solution from this study and the dashed curves were computed with RADMC-3D. The lower panels show the relative differences between the two codes. |
|
In the text |
Fig. 6 SED profiles for the axis-symmetric envelope with (blue curve) and (orange curve). The left and right panels correspond to i = 12.5 and 77.5 deg, respectively. The cross marks represent the solution from this study and the dashed curves were computed with RADMC-3D. The lower panels show the relative differences between the two codes. |
|
In the text |
Fig. 7 Images at λ = 2.3, 4.5, and 12.1 μm (left, middle, and right panels, respectively) of the 10 AU inner regions of the axis-symmetric envelope (Sect. 5.2) and computed from the DGFEM solution with the ray-tracing module (Appendix C). The top, middle and bottom panels correspond to the inclinations i = 12.5, 77.5, and 90 deg, respectively. The colour code shows the specific intensity value (in W m−2 Hz−1 sr−1) of the envelope , inside each pixel. The solid white lines show the iso-contours of . |
|
In the text |
Fig. 8 Image slices at x = 0 AU, for two inclinations i = 12.5 and 77.5 deg (left and right panels, respectively). The blue, orange, and green curves are the RADMC-3D intensities at λ = 2.3, 4.5, and 12.1 μm, respectively. The cross marks represent the solution from this study. The lower panels show the relative differences between the two codes. |
|
In the text |
Fig. 9 L2 norm of the relative differences in temperature between this study and RADMC-3D versus the number of elements in the computational grid (N = Nr × NΘ × Nμ × Nφ. The L2 norm was computed for the benchmark problem presented in Sec. 5.2. |
|
In the text |
Fig. 10 Emerging SEDs of RADMC-3D (solid black lines) and of the DGFEM code, computed from the post-processing ray-tracing procedure (orange dots), and from the DGFEM solution at the outer edge (blue crosses). The left panel corresponds to the optically thick spherically symmetric case (τ = 100, see Sect. 5.1) while the right panel corresponds to the optically thick disc benchmark from Sect. 5.2. The lower panels show the relative difference between the curves of this study and the benchmarks. |
|
In the text |
Fig. A.1 Geometry of a ray in the inner cavity. Point B is the opposite point of A, in the direction Ω. |
|
In the text |
Fig. B.1 Example of an element Di,j,k,l (dashed lines) using na = nc = nd = 3 and nb = 2 points. The left and right panels correspond to 2D slices in the (r, Θ) and (μ, φ) planes, respectively (we display here the local coordinate system given by Eq. B.2). The black dots correspond to the nodes where the solution is computed while the crosses correspond to the interpolated value of I used to compute the numerical flux at the interface, along the Θ coordinate. |
|
In the text |
Fig. C.1 Example of a ray (red line) normal to the image plane, crossing the spherical grid. In this example, we display a square image of size L × L and a circular image of radius R. The red crosses represent the intersections between the ray and the grid. The values of and Sv at these intersections are linearly interpolated from the grid adjacent values (red dots). |
|
In the text |
Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.