Abstract
We present an asymptotic preserving method for the radiative transfer equations in the framework of \(P_N\) method. An implicit and explicit numerical scheme is proposed to solve the \(P_N\) system based on the order analysis of the expansion coefficients of the specific intensity, where the order of each expansion coefficient is derived by the Chapman-Enskog method. The coefficients at higher-order are treated explicitly while those at lower-order are treated implicitly in each equation of the \(P_N\) system. Energy inequality is proved for this numerical scheme. Several numerical examples validate the efficiency of this scheme in both optically thick and thin regions.
Similar content being viewed by others
Data Availability
Enquiries about data availability should be directed to the authors.
References
Brunner, T. A.: Forms of approximate radiation transport. Sandia report, (2002)
Crestetto, A., Crouseilles, N., Dimarco, G., Lemou, M.: Asymptotically complexity diminishing schemes (ACDS) for kinetic equations in the diffusive scaling. J. Comput. Phys. 394, 243–262 (2019)
Densmore, J.: Asymptotic analysis of the spatial discretization of radiation absorption and re-emission in Implicit Monte Carlo. J. Comput. Phys. 230(4), 1116–1133 (2011)
Densmore, J., Park, H., Wollaber, A., Rauenzahn, R., Knoll, D.: Monte Carlo simulation methods in moment-based scale-bridging algorithms for thermal radiative-transfer problems. J. Comput. Phys. 284, 40–58 (2015)
Di, Y., Fan, Y., Kou, Z., Li, R., Wang, Y.: Filtered hyperbolic moment method for the Vlasov equation. J. Sci. Comput. 79(2), 969–991 (2019)
Discacciati, N., Hesthaven, J., Deep, R.: Controlling oscillations in high-order Discontinuous Galerkin schemes using artificial viscosity tuned by neural networks. J. Comput. Phys. 409, 109304 (2020)
Fleck, J., Cummings, J.: An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport. J. Comput. Phys. 8(3), 313–342 (1971)
Gentile, N.: Implicit Monte Carlo diffusion-an acceleration method for Monte Carlo time-dependent radiative transfer simulations. J. Comput. Phys. 172(2), 543–571 (2001)
Hammer, H., Park, H., Chacón, L.: A multi-dimensional, moment-accelerated deterministic particle method for time-dependent, multi-frequency thermal radiative transfer problems. J. Comput. Phys. 386, 653–674 (2019)
Hou, T., Li, R.: Computing nearly singular solutions using pseudo-spectral methods. J. Comput. Phys. 226(1), 379–397 (2007)
Jang, J., Li, F., Qiu, J., Xiong, T.: Analysis of asymptotic preserving dg-imex schemes for linear kinetic transport equations in a diffusive scaling. SIAM J. Numer. Anal. 52(4), 2048–2072 (2014)
Jang, J., Li, F., Qiu, J., Xiong, T.: High order asymptotic preserving DG-IMEX schemes for discrete-velocity kinetic equations in a diffusive scaling. J. Comput. Phys. 281, 199–224 (2015)
Jin, S., Levermore, C.: The discrete-ordinate method in diffusive regimes. Transp. Theory Stat. Phys 20(1–2), 413–439 (1991)
Jin, S., Levermore, C.: Fully discrete numerical transfer in diffusive regimes. Transp. Theory Stat. Phys 22(6), 739–791 (1993)
Jin, S., Pareschi, L., Toscani, G.: Uniformly accurate diffusive relaxation schemes for multiscale transport equations. SIAM J. Numer. Anal. 38(3), 913–936 (2000)
Kershaw, D.: Flux limiting nature’s own way. Technical Report UCRL-78378, Lawrence Livermore National Laboratory, Livermore, CA, (1976)
Klar, A.: An asymptotic-induced scheme for nonstationary transport equations in the diffusive limit. SIAM J. Numer. Anal 35(6), 1073–1094 (1998)
Koch, R., Krebs, W., Wittig, S., Viskanta, R.: Discrete ordinates quadrature schemes for multidimensional radiative transfer. J. Quant. Spectrosc. Ra. 53(4), 353–372 (1995)
Los Alamos National Laboratory. An implicit Monte Carlo code for thermal radiative transfer: Capabilities, development, and usag. LA-14195-MS, 2000
Laboure, V., McClarren, R., Hauck, C.: Implicit filtered \(P_N\) for high-energy density thermal radiation transport using discontinuous galerkin finite elements. J. Comput. Phys. 321, 624–643 (2016)
Laiu, M., Frank, M., Hauck, C.: A positive asymptotic-preserving scheme for linear kinetic transport equations. SIAM J. Sci. Comput. 41, A1500–A1526 (2019)
Larsen, A., Morel, J.: Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes. J. Comput. Phys. 69(2), 283–324 (1987)
Larsen, A., Morel, J.: Asymptotic solutions of numerical transport problems in optically thick, diffusive regimes ii. J. Comput. Phys. 83(1), 212–236 (1989)
Larsen, E., Kumar, A., Morel, J.: Properties of the implicitly time-differenced equations of thermal radiation transport. J. Comput. Phys. 238, 82–96 (2013)
Lathrop, K., Garlson, B.: Discrete ordinates angular quadrature of the neutron transport equation. Los Alamos Scientific Laboratory, (1965)
Lemou, E., Mieussens, L.: A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit. SIAM J. Sci. Comput. 31, 334–368 (2010)
Lewis, E., Miller, W.: Computational Methods in Neutron Transport. United States, (1993)
Li, W., Liu, C., Zhu, Y., Zhang, J., Xu, K.: Unified gas-kinetic wave-particle methods iii: Multiscale photon transport. J. Comput. Phys. 408, 109280 (2020)
Maginot, P., Ragusa, J., Morel, J.: High-order solution methods for grey discrete ordinates thermal radiative transfer. J. Comput. Phys. 327, 719–746 (2016)
Mathews, K.: On the propagation of rays in discrete ordinates. Nucl. Sci. Eng. 132, 155–180 (1999)
McClarren, R., Evans, T., Lowrie, R., Densmore, J.: Semi-implicit time integration for PN thermal radiative transfer. J. Comput. Phys. 227(16), 7561–7586 (2008)
McClarren, R., Hauck, C.: Robust and accurate filtered spherical harmonics expansions for radiative transfer. J. Comput. Phys. 229(16), 5597–5614 (2010)
McClarren, R., Hauck, C.: Simulating radiative transfer with filtered spherical harmonics. Phys. Lett. A 374(22), 2290–2296 (2010)
McClarren, R., Holloway, J., Brunner, T.: On solutions to the \(P_n\) equations for thermal radiative transfer. J. Comput. Phys. 227(5), 2864–2885 (2008)
Mieussens, L.: On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models. J. Comput. Phys. 253, 138–156 (2013)
Morel, J., Wareing, T., Lowrie, R., Parsons, D.: Analysis of ray-effect mitigation techniques. Nucl. Sci. Eng. 144, 1–22 (2003)
Olson, G.: Second-order time evolution of \(P_N\) equations for radiation transport. J. Comput. Phys. 228(8), 3072–3083 (2009)
Park, H., Knoll, D., Rauenzahn, R., Wollaber, A., Densmore, J.: A consistent, moment-based, multiscale solution approach for thermal radiative transfer problems. Transp. Theory Stat. Phys. 41(3–4), 284–303 (2012)
Peng, Z., Cheng, Y., Qiu, J., Li, F.: Stability-enhanced AP IMEX-LDG schemes for linear kinetic transport equations under a diffusive scaling. J. Comput. Phys. 415, 109485 (2020)
Seibold, B., Frank, M.: Starmap code. website. http://www.math.temple.edu/~seibold/research/starmap
Seibold, B., Frank, M.: Starmap-a second order staggered grid method for spherical harmonics moment equations of radiative transfer. ACM T. Math. Software (TOMS) 41(1), 4 (2014)
Shi, Y., Song, P., Sun, W.: An asymptotic preserving unified gas kinetic particle method for radiative transfer equations. J. Comput. Phys. 420, 109687 (2020)
Sun, W., Jiang, S., Xu, K.: An asymptotic preserving unified gas kinetic scheme for gray radiative transfer equations. J. Comput. Phys. 285(15), 265–279 (2015)
Sun, W., Jiang, S., Xu, K.: An asymptotic preserving implicit unified gas kinetic scheme for frequency-dependent radiative transfer equations. Int. J. Numer. Anal. Mod. 15(1–2), 134–153 (2018)
Tang, M., Wang, L., Zhang, X.: Accurate front capturing asymptotic preserving scheme for nonlinear gray radiative transfer equation. SIAM J. Sci. Comput. 43(3), B759–B783 (2021)
Warsa, J., Wareing, T., Morel, J.: Krylov iterative methods and the degraded effectiveness of diffusion synthetic acceleration for multidimensional \(S_N\) calculations in problems with material discontinuities. Nucl. Sci. Eng. 147, 218–248 (2004)
Xiong, T., Jang, J., Li, F., Qiu, J.: High order asymptotic preserving nodal discontinuous Galerkin IMEX schemes for the BGK equation. J. Comput. Phys. 284, 70–94 (2015)
Xiong, T., Sun, W., Shi, Y., Song, P.: High order asymptotic preserving discontinuous Galerkin methods for gray radiative transfer equations. arXiv:2011.14090, (2020)
Yan, B., Jin, S.: A successive penalty-based asymptotic-preserving scheme for kinetic equations. SIAM J. Sci. Comput. 35(1), A150–A172 (2013)
Acknowledgements
We thank Prof. Ruo Li from PKU, Prof. Zhenning Cai from NUS, Prof. Tao Xiong from XMU, Prof. Kailiang Wu from SUSTech, Dr. Zhichao Peng from MSU and Prof. Jiequan Li, Prof. Wenjun Sun, Dr. Yi Shi from IAPCM for their valuable suggestions.
Funding
Weiming Li is partially supported by the Science Challenge Project (No. TZ2016002) and the National Natural Science Foundation of China (12001051). Peng Song is partially supported by the Science Challenge Project (No. TZ2016002), the CAEP foundation (No.CX20200026). The work of Yanli Wang is partially supported by Science Challenge Project (No. TZ2016002) and the National Natural Science Foundation of China (Grant No. 12171026, U1930402 and 12031013).
Author information
Authors and Affiliations
Corresponding author
Ethics declarations
Conflict of interest
The authors have not disclosed any competing interest.
Additional information
Publisher's Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix
Appendix
1.1 The Gray Approximation of the Radiative Transfer Equations for 1D Angle Problem and Related \(P_N\) Equations
The time-dependent gray approximation of the radiative transfer equations [31] in a one-dimensional planar geometry medium has the form as
where \(I(t, x, \mu )\) is the specific intensity of radiation, \(\mu = \cos \theta \in [-1, 1]\) is the internal coordinate associated with the angle \(\theta \in [0, \pi ]\). T(t, x) is the material temperature, and \(\sigma \) is the absorption opacity. Moreover, the one-dimensional form of (2.10) is reduced into
For (A.1), the basis function for the \(P_N\) method is the Legendre polynomials. The moments are defined as
where \(P_l\) is the Legendre polynomial. Then, the \(P_N\) equations for (A.1) are
where \(\bar{\mathbf{I}}= (I_0, I_1, \cdots , I_M)\) and \(e_1 = (1, 0, \cdots , 0)^T\). \(\mathbf{B}^{\mathrm{low}}\) and \(\mathbf{B}^{\mathrm{up}}\) are triangular matrix with the non-zero entries as
1.2 Inflow Boundary Condition for PN Equations
We implement the inflow boundary condition for the \(P_N\) equations by specifying the values of coefficients of the \(P_N\) system in ghost cells. Choosing the left boundary as an example, the incoming specific intensity incident on the boundary interface is
For the \(P_N\) method, the numerical boundary can be rewritten as
where \(I^{i}(\mu )\) is the specific intensity at the left boundary of the area. Then, the expansion coefficient at the ghost cell is
The implementation of the inflow boundary condition in 2D is similar in spirit to that of 1D. Supposing \({\varvec{n}}\) is the outward normal of the boundary interface, the incident specific intensity on the boundary is
For the \(P_N\) method, the numerical boundary can be rewritten as
where \(I^{i}({\varvec{\Omega }})\) is the specific intensity on the interior side of the boundary interface. Thus, the expansion coefficient at the ghost cell is
The integration \(\int _{{\varvec{\Omega }}\cdot {\varvec{n}}> 0} Y^k_j({\varvec{\Omega }}) \overline{Y^m_l}({\varvec{\Omega }}) \,\mathrm {d}{\varvec{\Omega }}\) does not depend on the specific numerical solutions, and is pre-computed.
1.3 Proof of Proposition 2
In this section, the proof of the Proposition 2 is proposed here.
Proof
Following the method in [39], we will begin the Fourier analysis of (3.18) for the \(P_1\) system of the linear equation system (A.2). The result can be extended to the generalized \(P_N\) system naturally. We will first discuss two special cases where \(\xi = 0, \pi \). Therein, \(\mathbf{C}\) is reduced into a real diagonal matrix with the maximum eigenvalues equaling 1. According to the principle, the numerical scheme is stable.
Then, we study the general case by considering two scenarios according to the time step length:
-
1.
\(\epsilon < \Delta x\), where
$$\begin{aligned} \Delta t = C \Delta x^2. \end{aligned}$$(A.12)Substituting the time step length (A.12) into (3.36), we can find that \(\lambda _i, i=1,2\) are functions of C, \(\frac{\epsilon }{\Delta x}\), \(\alpha \) and \(\xi \). Introducing two variables as \(\beta _1 = \log _{10}(C)\) and \(\beta _2 = \log _{10} (\epsilon / \Delta x)\), the stability regions are plotted in Fig. 13 with fixed \(\alpha \). Here the discrete wave number \(\xi \) is uniformly taken from \([0, 2\pi ]\) with 200 samples. In this case, due to the definition of C which is the CFL number and \(\epsilon < \Delta x\), the range for \(\beta _1\) and \(\beta _2\) is changed into
$$\begin{aligned} \beta _1< 0, \qquad \beta _2 < 0. \end{aligned}$$(A.13)Moreover, it is natural to demand that \(\epsilon< \Delta x < 0.4\). Thus, \(\alpha \) is taken uniformly from \([0, \exp (-1/0.16)]\) with 100 samples and six cases are shown in Fig. 13 due to their similar behavior. From Fig. 13, we can find that when \(\alpha = 0\), the numerical scheme is always stable. However, with the increase of \(\alpha \), the stability region is becoming smaller, especially when the CFL number C is large and the radio \(\epsilon / \Delta x\) is small. We find that when \(\alpha = \exp (-1/0.16)\), the numerical scheme is stable when \(\log _{10}(\epsilon / \Delta x) > -2.5\). Noting that when \(\alpha = \exp (-1/0.16)\), which means \(\epsilon = 0.4\), \(\log _{10}(\epsilon / \Delta x)\) is always larger than \(-2.5\) for \(\Delta x < 1\). This indicates that in the simulation of benchmark problems, the stability condition is always satisfied.
-
2.
\(\epsilon > \Delta x\), where
$$\begin{aligned} \Delta t = C \epsilon \Delta x. \end{aligned}$$(A.14)Substituting the time step length (A.14) into (3.36), we can easily find that \(\lambda _i, i = 1, 2\) are also the function of C, \(\frac{\epsilon }{\Delta x}\), \(\alpha \) and \(\xi \). Introducing the same two variables \(\beta _{i}, i = 1,2\), we plot the stability regions in Fig. 14 with fixed \(\alpha \). Here the discrete wave number \(\xi \) is uniformly taken from \([0, 2\pi ]\) with 200 samples. Since
$$\begin{aligned} \alpha = \exp \left( -\frac{1}{\epsilon ^2}\right) , \end{aligned}$$(A.15)and assuming \(\epsilon < 1\) in the numerical test, \(\alpha \) is taken uniformly from [0, 0.5] with 100 samples. As their behavior is similar, the six cases \(\alpha = 0, 0.05, 0.1, 0. 2, 0.3\) and 0.5 are plotted here to illustrate the result. Moreover, due to the definition of C which is the CFL number, and the condition that \(\epsilon > \Delta x\), it holds that
$$\begin{aligned} \beta _1 < 0, \qquad \beta _2 > 0. \end{aligned}$$(A.16)From Fig. 14, we can find that the numerical scheme is stable under the time step length (A.14). In the numerical tests, the upper bound of \(\epsilon \) is \(\epsilon = 10^5 \Delta x\), which is large enough for the computational parameter.
\(\square \)
1.4 Proof of Theorem 2
In this section, the proof of Theorem 2 is proposed here.
Proof
We will take \(M=2\) as an example, and it could be extended to the general case naturally. Moreover, without loss of generality, we set \(a = c = C_v = \sigma = 1\) in the proof. When \(M=2\), (4.10) is reduced into
For (A.17a), multiplying it by \(I_{0,j}^{n+1}\), we can get that
For (A.17b) and (A.17c), shifting it backward one time step and multiplying \(3I_{1,j}^{n}\) and \(5I_{2, j}^{n}\) respectively, we can derive that
Summing (A.18) and (A.19) over j, then it holds that
where
Then we will begin from the approximation of \(A_1\) and \(A_2\). With some arrangement and the periodic boundary condition, \(A_1\) is changed into
With the estimation that
Similarly, it also holds that
Let
then together with (A.21), (A.22), (A.23), (A.24) and (A.25), (A.20) is reduced into
with
If it holds for \(\beta _6\) that
with the time step length (3.32), then we can derive the stability result (4.16). Precisely, with (A.17a) and (A.17d), we can derive that
Multiplying (A.29) with \((T_{j}^4)^{n+1}\) and summing over j, it holds with (A.26)
We derive the energy stability (4.16). The only point left is to prove (A.28), which we will be done in two cases:
-
1.
\(\epsilon > \Delta x\), in which case,
$$\begin{aligned} \Delta t= C \epsilon \Delta x. \end{aligned}$$(A.31)Substituting (A.31) into (A.27), we can deduce that
$$\begin{aligned} \beta _6 = \frac{9C^2\epsilon ^2}{10}\left( \frac{1}{1 - \alpha C }\right) + 6\alpha \epsilon ^2C(\alpha C - 1) - 3C \Delta x \epsilon . \end{aligned}$$(A.32)Thus if
$$\begin{aligned} 0< C < \min \left( \frac{\epsilon }{\alpha \Delta x},\frac{10 \Delta x}{3 \epsilon + 10 \Delta x \alpha } \right) , \end{aligned}$$(A.33)it holds that \(\beta _6 \le 0\).
For the coefficients \(\beta _i^2, i = 1,\cdots 5\), it requires that \(\beta _i^2 > 0\). Thus, from (A.25), it demands that
$$\begin{aligned} \frac{\epsilon \Delta x}{\Delta t} - \alpha > 0. \end{aligned}$$(A.34)Substituting (A.31) into (A.34), we can obtain that
$$\begin{aligned} C < \frac{1}{\alpha }. \end{aligned}$$(A.35)Thus, the constrain on C is changed into
$$\begin{aligned} 0< C < \min \left( \frac{1}{\alpha },\frac{10 \Delta x}{3 \epsilon + 10 \Delta x \alpha } \right) . \end{aligned}$$(A.36) -
2.
\(\epsilon < \Delta x\), in which case
$$\begin{aligned} \Delta t= C \Delta x^2. \end{aligned}$$(A.37)Substituting (A.37) into (A.27), we can deduce that
$$\begin{aligned} \beta _6 = \frac{9}{10}C^2\epsilon \Delta x^2\left( \frac{1}{\epsilon - \alpha C \Delta x}\right) + 6\alpha \Delta x C(\alpha C \Delta x - \epsilon ) - 3C \Delta x^2. \end{aligned}$$(A.38)Thus if
$$\begin{aligned} 0< C < \min \left( \frac{\epsilon }{\alpha \Delta x},\frac{10 \epsilon }{3 \epsilon + 10 \Delta x \alpha } \right) , \end{aligned}$$(A.39)it holds that \(\beta _6 \le 0\). Similarly, we can verify that the constrain \(\beta _i^2 > 0, i = 1, \cdots 5\) will not affect the condition (A.39), then the proof is finished.
For \(\epsilon < \Delta x\), it is always true that \(\alpha = \exp (-1/\epsilon ^2)\) is quite small, and (A.39) could be reduced into
$$\begin{aligned} 0< C < \frac{10}{3}. \end{aligned}$$(A.40)
\(\square \)
1.5 Analysis of the Higher-Order Scheme
From the test of the AP property for the numerical scheme, we found that even for the IMEX3 scheme with WENO reconstruction, the convergence order is only two. Analysis of the numerical scheme shows that when solving \(T^{n+1}\), the fourth-order polynomial equation of \(T^{n+1}\) is solved, where \((T^{n+1})^{4}\) is approximated as
instead of
where \(T_{i}\) is the cell average of cell i. Noting that
and
thus, it holds
Therefore, the convergence order of the whole numerical scheme is at most two.
Rights and permissions
About this article
Cite this article
Fu, J., Li, W., Song, P. et al. An Asymptotic-Preserving IMEX Method for Nonlinear Radiative Transfer Equation. J Sci Comput 92, 27 (2022). https://doi.org/10.1007/s10915-022-01870-3
Received:
Revised:
Accepted:
Published:
DOI: https://doi.org/10.1007/s10915-022-01870-3