A New Spherical Harmonics Scheme for Multi-Dimensional Radiation Transport I: Static Matter Configurations

Recent work by McClarren & Hauck [29] suggests that the filtered spherical harmonics method represents an efficient, robust, and accurate method for radiation transport, at least in the two-dimensional (2D) case. We extend their work to the three-dimensional (3D) case and find that all of the advantages of the filtering approach identified in 2D are present also in the 3D case. We reformulate the filter operation in a way that is independent of the timestep and of the spatial discretization. We also explore different second- and fourth-order filters and find that the second-order ones yield significantly better results. Overall, our findings suggest that the filtered spherical harmonics approach represents a very promising method for 3D radiation transport calculations.Comment: 29 pages, 13 figures. Version matching the one in Journal of Computational Physic

A class of Galerkin schemes for time-dependent radiative transfer

The numerical solution of time-dependent radiative transfer problems is challenging, both, due to the high dimension as well as the anisotropic structure of the underlying integro-partial differential equation. In this paper we propose a general framework for designing numerical methods for time-dependent radiative transfer based on a Galerkin discretization in space and angle combined with appropriate time stepping schemes. This allows us to systematically incorporate boundary conditions and to preserve basic properties like exponential stability and decay to equilibrium also on the discrete level. We present the basic a-priori error analysis and provide abstract error estimates that cover a wide class of methods. The starting point for our considerations is to rewrite the radiative transfer problem as a system of evolution equations which has a similar structure like first order hyperbolic systems in acoustics or electrodynamics. This analogy allows us to generalize the main arguments of the numerical analysis for such applications to the radiative transfer problem under investigation. We also discuss a particular discretization scheme based on a truncated spherical harmonic expansion in angle, a finite element discretization in space, and the implicit Euler method in time. The performance of the resulting mixed PN-finite element time stepping scheme is demonstrated by computational results

Numerical analysis of a spherical harmonic discontinuous Galerkin method for scaled radiative transfer equations with isotropic scattering

In highly diffusion regimes when the mean free path $\varepsilon$ tends to zero, the radiative transfer equation has an asymptotic behavior which is governed by a diffusion equation and the corresponding boundary condition. Generally, a numerical scheme for solving this problem has the truncation error containing an $\varepsilon^{-1}$ contribution, that leads to a nonuniform convergence for small $\varepsilon$. Such phenomenons require high resolutions of discretizations, which degrades the performance of the numerical scheme in the diffusion limit. In this paper, we first provide a--priori estimates for the scaled spherical harmonic ($P_N$) radiative transfer equation. Then we present an error analysis for the spherical harmonic discontinuous Galerkin (DG) method of the scaled radiative transfer equation showing that, under some mild assumptions, its solutions converge uniformly in $\varepsilon$ to the solution of the scaled radiative transfer equation. We further present an optimal convergence result for the DG method with the upwind flux on Cartesian grids. Error estimates of $\left(1+\mathcal{O}(\varepsilon)\right)h^{k+1}$ (where $h$ is the maximum element length) are obtained when tensor product polynomials of degree at most $k$ are used

Incorporating reflection boundary conditions in the Neumann series radiative transport equation: Application to photon propagation and reconstruction in diffuse optical imaging

We propose a formalism to incorporate boundary conditions in a Neumann-series-based radiative transport equation. The formalism accurately models the reflection of photons at the tissue-external medium interface using Fresnel’s equations. The formalism was used to develop a gradient descent-based image reconstruction technique. The proposed methods were implemented for 3D diffuse optical imaging. In computational studies, it was observed that the average root-mean-square error (RMSE) for the output images and the estimated absorption coefficients reduced by 38% and 84%, respectively, when the reflection boundary conditions were incorporated. These results demonstrate the importance of incorporating boundary conditions that model the reflection of photons at the tissue-external medium interface

An algorithm for computing the 2D structure of fast rotating stars

Stars may be understood as self-gravitating masses of a compressible fluid whose radiative cooling is compensated by nuclear reactions or gravitational contraction. The understanding of their time evolution requires the use of detailed models that account for a complex microphysics including that of opacities, equation of state and nuclear reactions. The present stellar models are essentially one-dimensional, namely spherically symmetric. However, the interpretation of recent data like the surface abundances of elements or the distribution of internal rotation have reached the limits of validity of one-dimensional models because of their very simplified representation of large-scale fluid flows. In this article, we describe the ESTER code, which is the first code able to compute in a consistent way a two-dimensional model of a fast rotating star including its large-scale flows. Compared to classical 1D stellar evolution codes, many numerical innovations have been introduced to deal with this complex problem. First, the spectral discretization based on spherical harmonics and Chebyshev polynomials is used to represent the 2D axisymmetric fields. A nonlinear mapping maps the spheroidal star and allows a smooth spectral representation of the fields. The properties of Picard and Newton iterations for solving the nonlinear partial differential equations of the problem are discussed. It turns out that the Picard scheme is efficient on the computation of the simple polytropic stars, but Newton algorithm is unsurpassed when stellar models include complex microphysics. Finally, we discuss the numerical efficiency of our solver of Newton iterations. This linear solver combines the iterative Conjugate Gradient Squared algorithm together with an LU-factorization serving as a preconditionner of the Jacobian matrix.Comment: 40 pages, 12 figures, accepted in J. Comput. Physic