Unified Gas-kinetic Wave-Particle Methods III: Multiscale Photon Transport

TL;DR

This paper introduces an advanced multiscale photon transport method that seamlessly transitions between particle tracking and diffusion regimes, significantly reducing computational costs and memory usage while maintaining accuracy across regimes.

Contribution

The paper extends the UGKWP method to photon transport, unifying wave and particle descriptions for efficient multiscale simulations without ray effects or high noise.

Findings

Requires less memory than $S_n$-type UGKS

Reduces statistical noise compared to IMC

Achieves several-order-of-magnitude computational savings

Abstract

In this paper, we extend the unified gas-kinetic wave-particle (UGKWP) method to the multiscale photon transport. In this method, the photon free streaming and scattering processes are treated in an un-splitting way. The duality descriptions, namely the simulation particle and distribution function, are utilized to describe the photon. By accurately recovering the governing equations of the unified gas-kinetic scheme (UGKS), the UGKWP preserves the multiscale dynamics of photon transport from optically thin to optically thick regime. In the optically thin regime, the UGKWP becomes a Monte Carlo type particle tracking method, while in the optically thick regime, the UGKWP becomes a diffusion equation solver. The local photon dynamics of the UGKWP, as well as the proportion of wave-described and particle-described photons are automatically adapted according to the numerical resolution and transport regime. Compared to the $S_n$ -type UGKS, the UGKWP requires less memory cost and does not suffer ray effect. Compared to the implicit Monte Carlo (IMC) method, the statistical noise of UGKWP is greatly reduced and computational efficiency is significantly improved in the optically thick regime. Several numerical examples covering all transport regimes from the optically thin to optically thick are computed to validate the accuracy and efficiency of the UGKWP method. In comparison to the $S_n$ -type UGKS and IMC method, the UGKWP method may have several-order-of-magnitude reduction in computational cost and memory requirement in solving some multsicale transport problems.

Tables

Table 1. Table 1: Comparison of computational time of the Monte Carlo method and UGKWP for inflow into homogeneous medium under different parameters.

$ϵ$	MC	UGKWP
$1$	$0.12$s	$0.14$s
${10}^{-2}$	$1.49$s	$1.35$s
${10}^{-4}$	$2587$s	$16$s

Table 2. Table 2: Comparison of the maximum number of particles used by the Monte Carlo method and UGKWP for inflow into homogeneous medium under different parameters.

$ϵ$	MC	UGKWP
$1$	$13053$	$15612$
${10}^{-2}$	$8054$	$6920$
${10}^{-4}$	$8343$	$0$

Equations

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mathbf{\Omega }\cdot \nabla I={\sigma }_s\left(\frac{1}{4\pi }{\int }_{{\mathbb{S}}^2}I\mathrm{d}\mathbf{\Omega }-I\right)-{\sigma }_{a}I+G,$$

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mathbf{\Omega }\cdot \nabla I=\sigma \left(\frac{1}{4\pi }{\int }_{{\mathbb{S}}^2}I\mathrm{d}\mathbf{\Omega }-I\right),$$

$$\frac{{ϵ}^2}{c}\frac{\partial I}{\partial t}+ϵ\mathbf{\Omega }\cdot \nabla I=\sigma \left(\frac{1}{4\pi }E-I\right),$$

$$I_m^n=\frac{1}{V_m}{\int }_{D_m}I(t^n,\mathbf{x},\mathbf{\Omega })\mathrm{d}\mathbf{x}$$

$$E_m^n=\frac{1}{V_m}{\int }_{D_m}E(t^n,\mathbf{x})\mathrm{d}\mathbf{x}$$

$$\dfrac{I^{n+1}_{m}-I^{n}_{m}}{\Delta t}+\tsum\slimits@\limits_{\boldsymbol{\tau}_{j}\in\partial D_{m}}\dfrac{\phi_{\boldsymbol{\tau}_{j}}}{V_{m}}=\dfrac{c\sigma}{\epsilon^{2}}\left(E^{n+1}_{m}-I^{n+1}_{m}\right),$$

$$\dfrac{E^{n+1}_{m}-E^{n}_{m}}{\Delta t}+\tsum\slimits@\limits_{\boldsymbol{\tau}_{j}\in\partial D_{m}}\dfrac{\Phi_{\boldsymbol{\tau}_{j}}}{V_{m}}=0,$$

$${\varphi }_{{\mathbf{\tau }}_j}=\frac{c}{ϵ\Delta t}{\int }_{t^n}^{t^{n+1}}{\int }_{{\mathbf{\tau }}_j}\left(\mathbf{\Omega }\cdot {\mathbf{n}}_j\right)I(t,\mathbf{x},\mathbf{\Omega })\mathrm{d}l\mathrm{d}t,$$

$${\Phi }_{j+\frac{1}{2}}={\int }_{{\mathbb{S}}^2}{\varphi }_{{\mathbf{\tau }}_j}(\mathbf{\Omega })\mathrm{d}\mathbf{\Omega }.$$

$$\begin{array}{rl}I(t,\mathbf{x},\mathbf{\Omega })=& {\mathrm{e}}^{-\frac{c\sigma (t-t^n)}{{ϵ}^2}}I\left(t^n,\mathbf{x}-\frac{c}{ϵ}\mathbf{\Omega }(t-t^n)\right)\\ & +{\int }_{t^n}^t{\mathrm{e}}^{-\frac{c\sigma (t-s)}{{ϵ}^2}}\times \frac{c\sigma }{{ϵ}^2}\frac{1}{4\pi }E\left(s,\mathbf{x}-\frac{c}{ϵ}\mathbf{\Omega }(t-s)\right)\mathrm{d}s,\end{array}$$

$$\frac{I_j^{n+1}-I_j^n}{\Delta t}+\frac{1}{\Delta x}\frac{c}{ϵ}\left(\mathbf{\Omega }\cdot \mathbf{n}\right)\left(\left(I_j^n{1}_{\mathbf{\Omega }\cdot \mathbf{n}>0}+I_{j+1}^n{1}_{\mathbf{\Omega }\cdot \mathbf{n}<0}\right)-\left(I_{j-1}^n{1}_{\mathbf{\Omega }\cdot \mathbf{n}>0}+I_j^n{1}_{\mathbf{\Omega }\cdot \mathbf{n}<0}\right)\right)=0.$$

$$\frac{E_j^{n+1}-E_j^n}{\Delta t}-\frac{1}{\Delta x}\left(\frac{c}{3{\sigma }_{j+\frac{1}{2}}}\frac{E_{j+1}^n-E_j^n}{\Delta x}-\frac{c}{3{\sigma }_{j-\frac{1}{2}}}\frac{E_j^n-E_{j-1}^n}{\Delta x}\right)=0.$$

$$\frac{{ϵ}^2}{c}\frac{\partial I}{\partial t}+ϵ\mathbf{\Omega }\cdot \nabla I=\sigma \left(\frac{1}{4\pi }E-I\right).$$

$$I(t,\mathbf{x},\mathbf{\Omega })={\int }_0^t{e}^{-\frac{c\sigma }{{ϵ}^2}(t-s)}\frac{c\sigma }{{ϵ}^2}\frac{1}{4\pi }E\left(s,\mathbf{x}-\frac{c}{ϵ}\mathbf{\Omega }(t-s),\mathbf{\Omega }\right)ds+e^{-\frac{c\sigma }{{ϵ}^2}t}{I}_0\left(\mathbf{x}-\frac{c}{ϵ}\mathbf{\Omega }t\right).$$

$$E\left(s,\mathbf{x}-\frac{c}{ϵ}\mathbf{\Omega }(t-s),\mathbf{\Omega }\right)=E(t,\mathbf{x},\mathbf{\Omega })+E_{\mathbf{x}}(t,\mathbf{x},\mathbf{\Omega })\frac{c}{ϵ}\mathbf{\Omega }(s-t)+E_t(t,\mathbf{x},\mathbf{\Omega })(s-t)+O(\Delta {\mathbf{x}}^2,\Delta t^2),$$

$$I(t,\mathbf{x},\mathbf{\Omega })=\left(1-e^{-\frac{c\sigma }{{ϵ}^2}t}\right)\frac{1}{4\pi }{E}^e(t,\mathbf{x},\mathbf{\Omega })+e^{-\frac{c\sigma }{{ϵ}^2}t}{I}_0(\mathbf{x}-\frac{c}{ϵ}\mathbf{\Omega }t).$$

$$E^e(t,\mathbf{x},\mathbf{\Omega })=E(t,\mathbf{x},\mathbf{\Omega }),$$

$$E^e(t,\mathbf{x},\mathbf{\Omega })=E(t,\mathbf{x},\mathbf{\Omega })+\frac{e^{-\frac{c\sigma }{{ϵ}^2}t}\left(t+\frac{{ϵ}^2}{c\sigma }\right)-\frac{c\sigma }{{ϵ}^2}}{1-e^{-\frac{c\sigma }{{ϵ}^2}t}}\left(E_t(t,\mathbf{x},\mathbf{\Omega })+\frac{c}{ϵ}\Omega E_{\mathbf{x}}(t,\mathbf{x},\mathbf{\Omega })\right).$$

$$t_c=-\frac{{ϵ}^2}{c\sigma }\ln \xi ,$$

$$t_c\ge \Delta t=t^{n+1}-t^n$$

$${\mathbf{x}}^{n+1}={\mathbf{x}}^n+\frac{c{\mathbf{\Omega }}^n}{ϵ}\Delta t.$$

$$E_{m}^{n+1}=E_{m}^{n}-\dfrac{\Delta t}{V_{m}}\tsum\slimits@\limits_{\boldsymbol{\tau}_{j}\in\partial D_{m}}\Phi_{E,\boldsymbol{\tau}_{j}}-\dfrac{1}{V_{m}}\tsum\slimits@\limits_{\boldsymbol{\tau}_{j}\in\partial D_{m}}W_{\boldsymbol{\tau}_{j}},$$

$${\Phi }_{E,{\mathbf{\tau }}_j}=\frac{c}{ϵ\Delta t}{\int }_{t^n}^{t^{n+1}}{\int }_{t^n}^t{\int }_{{\mathbf{\tau }}_j}{\int }_{{\mathbb{S}}^2}\left(\mathbf{\Omega }\cdot {\mathbf{n}}_j\right)\frac{c{\sigma }_{{\mathbf{\tau }}_j}}{{ϵ}^2}{\mathrm{e}}^{-\frac{c{\sigma }_{{\mathbf{\tau }}_j}(t-s)}{{ϵ}^2}}\frac{1}{4\pi }E\left(s,\mathbf{x}-\frac{c}{ϵ}(t-s)\mathbf{\Omega }\right)\mathrm{d}\mathbf{\Omega }\mathrm{d}l\mathrm{d}s\mathrm{d}t.$$

$$\begin{array}{rl}E\left(s,{\mathbf{x}}_{{\mathbf{\tau }}_j}-\frac{c}{ϵ}(t-s)\mathbf{\Omega }\right)=& E(t^n,{\mathbf{x}}_{{\mathbf{\tau }}_j})+\frac{\partial E}{\partial t}(t^n,{\mathbf{x}}_{{\mathbf{\tau }}_j})\times (s-t_n)\\ & -\frac{c}{ϵ}(t-s)\mathbf{\Omega }\cdot \nabla E(t,{\mathbf{x}}_{{\mathbf{\tau }}_j}),\end{array}$$

$$\nabla E(t^n,{\mathbf{x}}_{{\mathbf{\tau }}_j})\approx \frac{E_m^n-E_k^n}{{\mathbf{x}}_m-{\mathbf{x}}_k},$$

$${\Phi }_{E,{\mathbf{\tau }}_j}=l_j{\alpha }_{{\mathbf{\tau }}_j}{\mathbf{n}}_j\cdot \left(\frac{E_m^n-E_k^n}{{\mathbf{x}}_m-{\mathbf{x}}_k}\right),$$

$${\alpha }_{{\mathbf{\tau }}_j}=-\frac{c}{3{\sigma }_{{\mathbf{\tau }}_j}}\left(-\frac{2}{\left(\frac{c{\sigma }_{{\mathbf{\tau }}_j}\Delta t}{{ϵ}^2}\right)}\left(1-{\mathrm{e}}^{-\frac{c{\sigma }_{{\mathbf{\tau }}_j}\Delta t}{{ϵ}^2}}\right)+\left(1+{\mathrm{e}}^{-\frac{c{\sigma }_{{\mathbf{\tau }}_j}\Delta t}{{ϵ}^2}}\right)\right).$$

$$W_{{\mathbf{\tau }}_j}=\text{tsum}{\text{slimits@}}_{k\in S({\mathbf{\tau }}_j)}\text{sign}({\mathbf{n}}_j\cdot {\mathbf{\Omega }}_k){\omega }_k,$$

$${\Phi }_{f,{\mathbf{\tau }}_j}=\frac{c}{ϵ\Delta t}{\int }_{t^n}^{t^{n+1}}{\int }_{{\mathbf{\tau }}_j}{\int }_{{\mathbb{S}}^2}\left(\mathbf{\Omega }\cdot {\mathbf{n}}_j\right){\mathrm{e}}^{-\frac{c{\sigma }_{{\mathbf{\tau }}_j}(t-t^n)}{{ϵ}^2}}I\left(t^n,\mathbf{x}-\frac{c}{ϵ}(t-t^n)\mathbf{\Omega },\mathbf{\Omega }\right)\mathrm{d}\mathbf{\Omega }\mathrm{d}l\mathrm{d}t.$$

$$E^{n+1}_{m}=E^{n}_{m}-\dfrac{\Delta t}{V_{m}}\tsum\slimits@\limits_{\boldsymbol{\tau}_{j}\in\partial D_{m}}l_{j}\alpha_{\boldsymbol{\tau}_{j}}\mathbf{n}_{j}\cdot\left(\dfrac{E^{n}_{m}-E^{n}_{k}}{\mathbf{x}_{m}-\mathbf{x}_{k}}\right)-\dfrac{1}{V_{m}}\tsum\slimits@\limits_{\boldsymbol{\tau}_{j}\in\partial D_{m}}W_{\boldsymbol{\tau}_{j}}.$$

Full text

††journal: Journal of LaTeX Templates

1 Introduction

The radiative transfer equation describes photon propagation in the background medium and has important applications in the fields of astrophysics [1], atmospheric physics [2], optical imaging [3] and so on. In this paper, we focus on the gray radiative transfer equation with isotropic scattering, which reads

[TABLE]

where $I(t,\mathbf{x},\boldsymbol{\Omega})$ is the specific intensity in time $t$, space $\mathbf{x}\in\mathbb{R}^{3}$, and angle $\boldsymbol{\Omega}$. Here $c$ is the speed of light, $\sigma_{s}$ is the scattering coefficient, $\sigma_{a}$ is the absorption coefficient, and $G$ is an internal source of photons.

There are typically two categories of numerical methods for solving the radiative transfer equations. The first category consists of the deterministic methods with different ways of discretizing and modeling, such as the discrete ordinate method (DOM) [4, 5, 6, 7] and the moment methods [8, 9, 10, 11]. The second category is the stochastic methods, for example, the Monte Carlo method [12, 13, 14]. The Monte Carlo method is a very popular method for solving the radiative transfer problems. In comparison with the deterministic methods, it is more efficient in optically thin regime especially for the multidimensional cases, and it does not suffer from the ray effect. However, as a particle method, it unavoidably has statistical noise. Also, it becomes inefficient when it comes to diffusive regime. In diffusive regime where the mean free path is small, photons may go through a huge number of scatterings during their lifetimes. Direct simulation of each scattering process for all particles makes the Monte Carlo method very expensive for capturing the diffusive solutions. On the other hand, in the diffusive regime the photon transport process could be well described by the diffusion equation, which could be solved efficiently. Based on this observation, many hybrid methods have been developed in order to improve the overall efficiency in different regimes [15, 16, 17, 18], where the Monte Carlo method is used in the optically thin regions and the diffusion equation is applied to the optically thick regions. However, as far as we know, there is still no unified principle for accurate domain decomposition for different regimes.

Another approach towards releasing the stiffness issue in the diffusive regime is to develop asymptotic-preserving (AP) DOM [19, 20, 21, 22, 23, 24, 25, 26, 27, 28] . One of the examples is the unified gas-kinetic scheme (UGKS), which couples the particles’ transport and collision process using a multiscale flux function obtained from the integral solution of the kinetic model equation. The cell size and time step are not restricted by the mean free path and mean collision time. It was developed initially in the field of rarefied gas dynamics [29, 30, 31] and has been applied to the field of radiative transfer [23, 24, 25, 26, 27, 28], plasma transport [32], and disperse multi-phase flow [33]. Since it is a discrete ordinate based numerical scheme, it has no statistical noise, but unavoidably suffers from the ray effect.

The UGKS provides a general methodology to construct multiscale simulation methods for the transport equations [30, 34]. It consists of two governing equations for the microscopic distribution function and macroscopic flow variables on the mesh size and time step scales, and a multiscale evolution solution for the interface distribution function. The time-dependent evolution solution at a cell interface covers the dynamics from the particle free transport to the hydrodynamic wave interaction with the variation of the ratio of the time step to the local particle collision time. The original UGKS is constructed based on the discrete velocity method (DVM) or DOM formulations [29]. However, the corresponding purely particle version and wave-particle version of UGKS can be constructed as well [35, 36]. In this work, we develop a novel unified gas-kinetic wave-particle (UGKWP) method for the multiscale photon transport, which combines the advantages of UGKS and the particle method. To facilitate understanding, we first describe a purely particle-based unified gas-kinetic particle (UGKP) method. In the UGKP method, the photons are described by the particle transport and collision, and this process is controlled by a multiscale integral formulation. More specifically, the Monte Carlo particle model is used to discretize the angular direction of the photon’s movement. Based on the particles’ transport nature in the discretized physical space, particles are categorized into two groups. Given a fixed time step, the free-stream particles are accurately tracked by following the trajectories of the simulation particles, while those particles that get scattered within the given time step are eliminated and re-sampled according to their macroscopic variables at the new time level. The fluxes across a cell interface from different type particles are taken into account to update cell averaged macroscopic variables. In such a way, the multiscale process is preserved by coupling particle free streaming and collision. Based on UGKP, a more efficient UGKWP method is proposed. Instead of representing all photons by simulation particles, a proportion of photon is represented by an analytical distribution function in UGKWP. Therefore, part of the evaluation of the fluxes, which were computed by particles in the UGKP method, could be done analytically with significant reduction in computational cost and statistical noise, especially in the diffusive regime. The multiscale flux function of the UGKS is precisely preserved in the UGKWP implementation. In the diffusive regime, the resulting algorithm would become a standard central difference scheme for the diffusion equation without any particle and DOM coordinate in the velocity space. In the optically thin regime, it gives a particle tracking method same as the Monte Carlo method. In the transition regime, the ratio of the time step over particle collision time determines the transport dynamics between the above two limits. The UGKWP is not a hybrid method of domain decomposition, but is a method to recover the multiscale modeling in UGKS [30] with the help of wave and particle in order to achieve a highest efficiency in the transport simulation in different regimes.

The rest of this paper is organized as follows. Section 2 briefly recalls the basic idea of the unified gas-kinetic scheme (UGKS) for the linear transport equation. Section 3 presents the UGKP method for linear photon transport. The UGKWP method is described in Section 4, while the extension to radiation-material coupling system is discussed in Section 5. In Section 6, numerical tests are presented to demonstrate the accuracy and efficiency of the UGKWP method. The final section is the conclusion.

2 Review of the UGKS for the linear transport equation

The unified gas-kinetic scheme (UGKS) was initially developed for the problems in the field of rarefied gas dynamics [29, 30], and has also been successfully applied to problems in radiative transfer under the finite volume framework [23, 24, 25, 26, 27, 28]. In this section, we review the basic idea of the UGKS using the example of the multidimensional linear transport equation in a purely scattering medium.

Consider

[TABLE]

which gives a non-dimensional equation

[TABLE]

where $E=\int_{\mathbb{S}^{2}}I(\boldsymbol{\Omega})\,{\rm d}\boldsymbol{\Omega}$. The non-dimensionalization is used same as that in [23], and $\epsilon$ is the ratio between the typical mean free path and the macroscopic length scale.

The UGKS is based on a finite volume framework. Consider a single cell $m$. Denote $D_{m}$ as the area covered by this cell, $\boldsymbol{\tau}_{j}$ as the interfaces of its boundary $\partial D_{m}$, $\mathbf{n}_{j}$ as the outward normal of $\boldsymbol{\tau}_{j}$, and $V_{m}$ as the volume of cell $m$. Define

[TABLE]

to be the averaged specific intensity $I$ over the spatial cell, and

[TABLE]

to be the averaged energy density function $E$ over a spatial cell. Under the finite volume framework, the discretizations of the two fundamental governing equations for microscopic and macroscopic variables in UGKS are,

[TABLE]

and

[TABLE]

where the microscopic and macroscopic flux terms are respectively

[TABLE]

and

[TABLE]

Based on the above microscopic and macroscopic governing equations, the key ingredient of the UGKS is the construction of the multiscale flux function by adopting the integral solution of the kinetic model equation (3). The integral solution of equation (3) along the characteristic line gives

[TABLE]

which is used to construct the numerical fluxes in equation (4) and (5). The integral solution couples transport with particle collisions, and bridges the kinetic and the hydrodynamic scale dynamics. Note that the dynamics in the evolution equations of UGKS through equations (4), (5), and (8), is related to time step $\Delta t$. The ratio of the time step to the local particle collision time determines the transport regime of the UGKS.

The numerical flux for microscopic and macroscopic variable updates is based on the piecewise linear initial reconstruction of $I$ and $E$ at the beginning of each time step. The discretization of $I$ is based on the DOM in the original UGKS. It has been proved in [23] that when $\sigma$ equals [math], the UGKS tends to a finite volume scheme which is consistent with free transport solution. More specifically, for 1D case with uniform grid, the scheme becomes

[TABLE]

In the diffusion limit, with a uniform mesh the UGKS scheme becomes a standard central difference method for the limit diffusion equation as $\epsilon$ tends to [math]. Specifically, for 1D case with uniform grid, in the diffusive limit the scheme is

[TABLE]

More details on the asymptotic analysis of the UGKS for the radiative transfer equation can be found in [23] and [24].

Following the methodology of the UGKS, we will construct two particle-based algorithms with multiscale transport property for recovering transport physics from the kinetic scale to the hydrodynamic scale, i.e., the unified gas-kinetic particle (UGKP) method and the unified gas-kinetic wave-particle (UGKWP) method. The UGKP is a purely particle based method whose computational cost and statistical noise remain the same in all flow regimes. The UGKWP method is an improvement of the UGKP method such that the distribution function is represented by the combination of simulation particles and an analytical distribution function, which makes it more efficient and less noisy than the UGKP in the diffusive regime. In both algorithms, for the kinetic scale particle free transport, the particle trajectories are tracked precisely; for those particles suffering collisions, their macroscopic variables will be updated and used to re-sample these particles subsequently. The multiscale particle methods for equation (8) are constructed through the tracking and re-sampling of particles with the help of updated macroscopic variables. We will first introduce the UGKP method in the next section.

3 The unified gas-kinetic particle (UGKP) method

In this section, we present a unified gas-kinetic particle (UGKP) method. Following the direct modeling methodology of UGKS, the evolution of microscopic simulation particle is coupled with the evolution of macroscopic energy in UGKP method. The method presented in this and the next section is based on a single linear transport equation,

[TABLE]

In the following two subsections, we will introduce the multiscale evolution equations of particle and macroscopic energy. For simplicity, the method will be presented for the two-dimensional case on uniform Cartesian mesh. Its extension to non-uniform mesh and 3D is straightforward.

3.1 Multiscale Particle evolution

The evolution of the simulation particle follows the integral solution of the linear transport equation (9),

[TABLE]

The local equilibrium $E$ can be expanded as

[TABLE]

and the integral solution can be written as

[TABLE]

The first order approximation gives

[TABLE]

and the second order approximation gives

[TABLE]

The reformulated integral solution Eq.(12) is the multiscale governing equation for particle evolution in UGKP method, which implies that the simulation particle has a probability of $e^{-\frac{c\sigma}{\epsilon^{2}}t}$ to free stream and has a probability of $1-e^{-\frac{c\sigma}{\epsilon^{2}}t}$ to get scattered. The scattered particle follows a velocity distribution $E^{e}(t,x,\boldsymbol{\Omega})$, and the first order approximation Eq.(13) is used in this paper. Different from the traditional Monte Carlo method which tracks the trajectory of each particle as shown in Fig. 1, UGKP tracks particle trajectory until its first scattering time, and then evolves the velocity distribution function. The scattered particles will be temporarily removed and then get re-sampled at $t^{n+1}$ from their evolved distribution, as shown in Fig. 2. The free-streaming particles are named as the $E_{f}$ particles, and the scattered particles are named as the $E^{+}$ particles. The total energy of the $E^{+}$ particles is called $E^{+}$ energy. In UGKP we do not need to resolve each scattering process by dividing every time step into several sub time steps and therefore achieve a high efficiency compared to the traditional Monte Carlo method.

The first scattering time $t_{c}$ is defined as the time when the first scattering happens to the particle. According to Eq.(12), the time $t_{c}$ can be obtained for each particle as

[TABLE]

where $\xi$ is a random number which follows a uniform distribution on $(0,1)$. Simulation particles can be categorized into two groups based on $t_{c}$. The first group of simulation particle satisfies

[TABLE]

and does not suffer any collision during a time step $\Delta t$, and this particle position is updated by

[TABLE]

The second group of particle satisfies $t_{c}<\Delta t$, those particles suffer one or multiple times of scattering during $[t_{c},t^{n+1}]$, with the assumption $t^{n}=0$ for simplification. Therefore, these particles will only be accurately tracked for $t\in[t^{n},t_{c})$ and be re-sampled from $E^{+}(t^{n+1},\mathbf{x},\boldsymbol{\Omega})$ at $t^{n+1}$. $E^{+}(t^{n+1},\mathbf{x},\boldsymbol{\Omega})$ is obtained from the evolved macroscopic energy. In the UGKP method, the evolution of particle and macroscopic energy are closely coupled. The multiscale evolution equation of macroscopic energy is given in the next subsection.

3.2 Multiscale macroscopic energy evolution

The update of the macroscopic variable follows the finite volume framework

[TABLE]

where the flux $\Phi_{E}$ is the diffusive flux and $W$ is the contribution of the free streaming flux. Based on (12), the diffusive flux is computed same as that in the UGKS,

[TABLE]

We adopt the piecewise linear reconstruction for $E$,

[TABLE]

and the explicit central difference discretization of $\nabla E$,

[TABLE]

where $\mathbf{x}_{m}$, $\mathbf{x}_{k}$ are the barycenters of cells $m$ and $k$ with the same interface $\boldsymbol{\tau}_{j}$. Then the diffusive flux is given by

[TABLE]

where

[TABLE]

The free transport flux $W$ is contributed by the particle free streaming, which is computed by

[TABLE]

where $\omega_{k}$ and $\boldsymbol{\Omega}_{k}$ are respectively the weight and velocity of the particle $k$, and $S(\boldsymbol{\tau}_{j})$ is the set of all particles which are transported across interface $\boldsymbol{\tau}_{j}$ during the free streaming process. It recovers the other part of UGKS flux

[TABLE]

With the diffusive flux Eq. (17) and the free transport flux Eq. (18), the macroscopic energy $E^{n+1}$ can be updated by

[TABLE]

With the updated $E^{n+1}$, the energy $E^{+}$ of the scattered particles can be obtained and be used to re-generate particles,

[TABLE]

where $Q(m)$ is the set of $E_{f}$ particles in cell $m$ at $t^{n+1}$. In summary, the evolution of particle provides the free-stream flux for macroscopic energy evolution, and the updated energy $E^{+}$ provides the distribution for the corresponding particles. The evolution of microscopic and macroscopic quantities are closely coupled, and the algorithm of UGKP is shown in Fig. 4.

4 Unified gas-kinetic wave-particle (UGKWP) method

4.1 The multiscale evolution of particles and macroscopic

variable

In UGKP, the $E^{+}$ particles are re-sampled from an equilibrium distribution at $t=t^{n+1}$. According to the integral solution (10), $E^{+}$ particles have a probability $1-\mathrm{e}^{-\frac{c\sigma\Delta t}{\epsilon^{2}}}$ to get scattered in the next time step from $t^{n+1}$ to $t^{n+2}$. We define the scattered particles in $t\in[t^{n+1},t^{n+2}]$ as $E^{+}_{W}$ particles and the non-scattering particles as $E^{+}_{P}$ particles. As the $E^{+}_{W}$ particles will be deleted after $t_{c}$ within $t\in[t^{n+1},t^{n+2}]$, and their contribution to the free transport flux could be computed analytically, only the $E^{+}_{P}$ particles need to be re-sampled and they keep on free transport in the step from $t^{n+1}$ to $t^{n+2}$. For the $E^{+}_{W}$ particles, we only need to store the distribution function instead of re-sampling particles from it. Based on this observation, we propose a more efficient unified gas-kinetic wave-particle (UGKWP) method.

The unified gas-kinetic wave-particle(UGKWP) method improves UGKP in two aspects:

Only $E^{+}_{P}$ particles are sampled as shown in Fig. 5; 2. 2.

The free transport flux contributed by $E^{+}_{W}$ can be calculated analytically.

The finite volume energy evolution follows

[TABLE]

where the diffusive flux $\Phi_{E,\boldsymbol{\tau}_{j}}$ is computed using Eq. (17), and the free transport flux contributed by $E^{+}_{W}$ is calculated by

[TABLE]

The free transport flux $W_{\boldsymbol{\tau}_{j}}$ contributed from the particles in Eq. (22) is computed by particle tracking,

[TABLE]

where $S(\boldsymbol{\tau}_{j})$ is the set of all particles which are transported across the cell interface $\boldsymbol{\tau}_{j}$ during the free streaming process. In the particle free transport process, the re-sampled $E^{+}_{P}$ particles will always have $t_{c}>\Delta t$ and keep free transport without collision in the whole time step. The other particles from $E_{f}$ will transport according to their individual generated $t_{c}$. With the solution of $E^{n+1}$, Eq.(21) is used to get $E^{+}$, which is divided into

[TABLE]

and

[TABLE]

and only particles from $E^{+}_{P}$ will be generated and keep on free transport in the whole next time step.

The algorithm of UGKWP is shown in Fig. 6.

4.2 Properties of the UGKWP algorithm

The unified gas-kinetic wave-particle method satisfies the following properties:

The microscopic and macroscopic evolutions are consistently evaluated based on the wave-particle decomposition. The macroscopic energy density $E$ is the sum of the energy from the un-scattering particles and the scattering one $E^{+}$. The evolution of $E^{+}_{W}$ particles can be evaluated analytically. 2. 2.

In the diffusive limit, $\epsilon\rightarrow 0$ and $\mathrm{e}^{-\frac{\sigma_{\boldsymbol{\tau}_{j}}\Delta t}{\epsilon^{2}}}\rightarrow 0$, essentially no particle can be freely transported and survived within a time step. Therefore, no particle is re-sampled from the $E^{+}$ energy. The algorithm becomes a time-explicit central difference solver of the diffusion equation. For instance, in 1D case with uniform grid, the scheme tends to the following limiting equation

[TABLE]

Note that in the diffusion regime, different from DOM method there is no discrete velocity space discretization in UGKWP. 3. 3.

In the free transport limit, $\sigma\rightarrow 0$ and $\mathrm{e}^{-\frac{\sigma_{\boldsymbol{\tau}_{j}}\Delta t}{\epsilon^{2}}}\rightarrow 1$, each particle is traced exactly by free transport with probability 1. In this case, the UGKWP could recover the exact solution for individual particle.

5 Extension to the coupled system of gray radiative transfer and material energy evolution

This section extends UGKWP to solve the coupled system of gray radiative transfer and material temperature equation,

[TABLE]

Define $u_{r}=aT^{4}$ and $\beta=\dfrac{\partial u_{r}}{\partial T}$, then the second equation can be re-written as

[TABLE]

The implicit Monte Carlo method proposed by Fleck and Cummings in [12] has been shown to be an effective technique for solving non-linear, time-dependent, radiative transfer problems and is widely used in the radiative transfer community. Fleck’s implicit Monte Carlo method uses an effective scattering process to approximate the absorption and emission of radiation by the background medium. This treatment allows it to take larger time steps than that in a purely explicit method. Here the similar semi-implicit discretization for material temperature evolution will be employed. Specifically, Eq. (25) is discretized by

[TABLE]

which gives

[TABLE]

With the definition

[TABLE]

substituting Eq. (26) into Eq. (24) yields

[TABLE]

An operator splitting approach is used to solve the above system subsequently by the linear kinetic equation

[TABLE]

the radiation energy exchange,

[TABLE]

and the update of material energy. Here Eq. (27) is solved using the algorithm introduced in Section 4. The exact solution of Eq. (28) is

[TABLE]

where $I$ is solved from Eq. (27) by UGKWP. The exponential decay term in Eq. (29) is implemented by modifying the weight of particles, while the source term in Eq. (29) is added to $E^{+}$ energy. Afterwards, the energy change of particles is summed up, and the material temperature and $u_{r}$ are updated by energy conservation. The above UGKWP can be further improved if the technique in [25] is incorporated into the current scheme, where both macroscopic equations for the radiation energy and material energy are solved iteratively first before updating the radiation intensity $I$.

6 Numerical Experiments

In this section, we present six numerical examples to validate the proposed UGKWP method. In all these examples, the UGKWP method obtains results which are in good agreement with that of the Monte Carlo method. At the same time, the UGKWP is more efficient and less noisy compared with the Monte Carlo method in the near diffusive regime. All computations are performed in sequential code on a computer with Intel i7-8700K CPU and $64$ GB memory. In all simulations, the time step is determined by $\Delta t=\text{ CFL}\times\epsilon\Delta x/c$, with $\text{CFL}=0.4$ for 1D examples and $\text{CFL}=0.2$ for 2D examples.

6.1 Inflow into purely scattering homogeneous medium

We first consider the behaviour of the UGKWP method for purely scattering homogeneous medium. Tests in this section are for non-dimensional linear equation

[TABLE]

defined on the semi-infinite spatial domain $x\in[0,\infty)$ with an isotropic inflow condition imposed on the left boundary. For the numerical simulation, the spatial domain is taken to be $[0,1]$. Inflow boundary condition is imposed at $x=0$ with the incoming specific intensity $I(t,0,\mu)=\frac{1}{2}$. The initial value is $I(\mu)=0$ for all $x$.

The results of both the UGKWP method and the Monte Carlo method are obtained using $200$ grids in space. As we are targeting to develop a method that automatically bridges the optically thin and optically thick regimes, the parameters cover the rarefied ($\epsilon\gg\Delta x$), the intermediate ($\epsilon\approx\Delta x$), and the diffusive ($\epsilon\ll\Delta x$) regimes, as defined in [37]. When $\epsilon$ is small, i.e. in the diffusive regime, the current method can use a much larger cell size and time step than the particle mean free path and collision time. Also, due to the exponential factor in re-sampling, fewer particles are used in UGKWP than the Monte Carlo method when $\epsilon$ is small, and UGKWP becomes much more efficient and less noisy than the Monte Carlo method in those cases. Table 1 compares the computational time of the Monte Carlo method and the UGKWP method for different $\epsilon$, while Table 2 compares the maximum number of particles used. Notice that no particle is generated in UGKWP method when $\epsilon=10^{-4}$ and the method only solves a diffusion equation. For $\epsilon=10^{-4}$, the Monte Carlo method takes more than half an hour while the UGKWP method takes only $16$ seconds, showing that the UGKWP method is much more efficient than the Monte Carlo method in the diffusive regimes.

In Figure 7 the numerical results of the UGKWP method are compared with the solutions of the Monte Carlo method for different $\epsilon$.

In this example, the UGKWP method gives solutions which are almost identical with the Monte Carlo solution for all flow regimes. Also, note that the UGKWP solution is smooth for $\epsilon=10^{-4}$, while the Monte Carlo solution has statistical noise.

6.2 Inflow into purely scattering heterogeneous medium

We next consider the UGKWP method for purely scattering heterogeneous medium. For this example, the set-up for the computational domain as well as that of the initial and boundary conditions are the same as in the previous section. Tests in this section are for the dimensional linear equation

[TABLE]

Both the UGKWP method and the Monte Carlo method use $200$ grids in space.

We calculate two test cases. For the first test case, we take $\sigma=10000\text{arctan}(1-x)$. In this case, radiation passes through an optically thick medium near the left boundary and gets into the region with a gradually reducing of optically thickness. In Figure 8(a) we compare the solutions of the UGKWP method and the Monte Carlo method for $t=1000$ and they agree with each other closely. Also, in this case where the medium is relatively optically thick overall, the UGKWP produces solutions much smoother than those given by the Monte Carlo method. For this test case, the Monte Carlo method uses a maximum of $17972$ particles and takes $5127$ seconds. The UGKWP method uses a maxium of $85$ particles and takes $18$ seconds. Therefore, the UGKWP method is much more efficient than the Monte Carlo method for this test case.

For the second test case, we take $\sigma=100/\text{arcsin}(1-x)$, implying that the medium gets more and more optically thick from the left to right boundaries. In Figure 8(b) the numerical result of the UGKWP method is compared with the solution of the Monte Carlo method at $t=20$ and they are almost identical except for the statistical noise. The Monte Carlo method uses a maximum of $24598$ particles and spends $12$ seconds while the UGKWP method uses a maximum of $24823$ particles and spends $10$ seconds in computing this case. In this test case the medium is relatively optically thin overall, and the UGKWP method has similar time cost as the Monte Carlo method.

6.3 Marshak wave problem

This section studies the Marshak wave problem where radiation is coupled with material medium. Consider the system

[TABLE]

for semi-infinite domain $x\in[0,+\infty)$ with a computational one $x\in[0,0.5]$. We take absorption coefficient to be $\sigma=\dfrac{30}{T^{3}}$, the speed of light to be $c=29.98$, the parameter $a=0.01372$ and the specific heat to be $C_{v}=0.3$. The initial material temperature is set to be $10^{-2}$ and initially material and radiation energy are at equilibrium. A constant isotropic inflow radiation intensity with an equivalent radiation temperature of $1$ is imposed on the left boundary. Both the UGKWP method and the Monte Carlo method use $200$ grids in space. For this example, the implicit Monte Carlo method uses a maximum of $45277$ particles, while the UGKWP method uses a maximum of $42892$ particles. It takes the implicit Monte Carlo method $66$ seconds to compute until $t=1.0$ and the UGKWP method $60$ seconds.

In Figure 9(a), the computed radiation wave front at time $t=0.33$, $0.66$ and $1.0$ are given, while Figure 9(b) presents the computed material temperature. The solutions of the UGKWP method are shown to be consistent with those of the Monte Carlo method for both radiation and material temperature.

6.4 Line-source problem in purely scattering homogeneous

medium

The next problem we look at is the line-source problem in purely scattering medium. In this test case, we consider equation (3) for $\epsilon=1$ and $c=1$. The spatial 2D computational domain is $[-1.5,1.5]\times[-1.5,1.5]$. The initial density distribution is

[TABLE]

which means that initially all particles are concentrated at $x=y=0$ and they spread out over the time. This is a particularly difficult problem for the $S_{n}$ based methods because they suffer from severe ray effect. Previous studies have used this problem for comparisons between methods [38] and as a test case for schemes to mitigate the ray effect [39].

Both the UGKWP method and the Monte Carlo method use $201\times 201$ cells in space. The UGKWP method uses $1.633\times 10^{7}$ particles while the Monte Carlo method uses $1.638\times 10^{7}$ particles and the same spatial mesh as the UGKWP method. Because this example considers the kinetic regime, the number of particles used in the UGKWP method is almost the same as that used in the Monte Carlo method and the computation time is similar, with $1530$ seconds for the UGKWP method and $1034$ seconds for the Monte Carlo method. The UGKWP method is slightly more expensive due to its additional overheads.

Fig.10 compares the solution of $E$ between the UGKWP method, the Monte Carlo method and the $S_{8}$ method for time $t=1.2$. For the $S_{8}$ computation, we use an equal-weight quadrature set, $201\times 201$ cells for the spatial mesh and the DOM-UGKS. The UGKWP solution is consistent with the Monte Carlo solution. Both solutions display a sharp spike near the wave front and has a smooth nonzero region behind the wave front, and neither suffers from the ray effect. On the other hand, the $S_{8}$ solution is qualitatively incorrect due to the ray effect. This example demonstrates that the UGKWP method preserves the advantage of the Monte Carlo method in the rarefied regime without suffering from the ray effect.

6.5 Emitting isotropic source in lattice medium

In this section we study a problem with multiple medium [40] by considering the following equation

[TABLE]

The diffusion term $b\Delta T$ is solved by using a backward Euler method. The two-dimensional physical space $[0,7]\times[0,7]$ consists of a set of squares belonging to a strongly absorbing medium and a background of weakly scattering medium. The specific layout of the problem is given in Fig.11, where the light green regions and the white region are purely scattering with $\sigma_{s}=2$ and $\sigma_{a}=0$; the red regions are pure absorbers with $\sigma_{s}=0$ and $\sigma_{a}=2000$. In the white region in the center there is an isotropic source $G=\dfrac{1}{4\pi}$. The initial material temperature is $10^{-2}$ and initially the radiation and material temperature are at equilibrium. We take $C_{v}=0.3$, $b=0.03$ and $c=a=\epsilon=1$ to simulate the case.

Both the UGKWP method and the Monte Carlo method use a mesh of $280\times 280$ in physical space. The results at time $t=2$ are summarized in Fig. 12 with the natural logarithm of density distribution shown in contour and slice, and they are in agreement with each other. This is a test case which includes material in both the kinetic and the diffusive regimes, but most of the time particles travel in the kinetic regime. The Monte Carlo method uses $10^{8}$ particles and took $5692$ seconds, while the UGKWP method uses $8.5\times 10^{7}$ particles and took $3092$ seconds, showing the UGKWP method to be more efficient than the Monte Carlo method for this case where there is mixed medium.

6.6 A hohlraum problem

This examples considers a hohlraum problem similar to the one studied in [40]. We study equation (24) with $G=0$. The problem layout is given in Fig.13. This is a problem of purely absorbing medium where the absorption coefficient relies on the material temperature. As the material temperature varies, the values of the absorption coefficient cover a wide range, presenting challenges to numerical methods.

On the left boundary there is an isotropic source with the specific intensity $I(\boldsymbol{\Omega})=\dfrac{1}{4\pi}$. Initially radiation and material temperature are at equilibrium and the initial material temperature is $T=10^{-2}$.

The UGKWP method uses $400\times 400$ grids in space while the Monte Carlo method use $800\times 800$ grid in space. We compute for $t=2$. The Monte Carlo method uses a maximum of $2\times 10^{7}$ particles and takes $9.8$ hours. The UGKWP method uses a maximum of $3\times 10^{6}$ particles and takes $1406$ seconds, therefore the UGKWP method is $25$ times more efficient than the Monte Carlo method for this test case. For this test case, the initial absorption coefficient on the boundary is extremely large, and the Monte Carlo method took an especially long time over the first time step when the inflow radiation first heats up the boundary material. The UGKWP method, on the other hand, has no such problems.

The authors of [40] pointed out solution to this problem should contain two characteristics: the central block should be heated non-uniformly; and the region behind the block with respect to the source should have less radiation energy then regions within the line of sight of the source. Their studies showed that the diffusion equation could not capture these two characteristics. In Fig. 14 and 15, we compare the solution of the UGKWP method and the Monte Carlo method for energy density and material temperature at $t=2$ and they are all consistent with each other. Both solutions agree with each other such that there is less radiation to the right of the block. This shows that the UGKWP method keeps the accuracy of the Monte Carlo method and is more accurate than the diffusion approximation in problems which are essentially multiscale in nature.

7 Conclusion

In this paper, for the first time a unified gas-kinetic wave-particle method is proposed to simulate photon transport. Due to the particle wave decomposition and their dynamic coupling, the UGKWP method becomes a multiscale method to simulate transport physics with the most efficient way in different regimes. For the linear transport equation, this method recovers the solution of the diffusion equation in the optically thick limit without constraint on the time step being less than the photon’s mean collision time. At the same time, it gives the exact solution in the free transport regime. In the transition regime, both macroscopic wave and kinetic particle contribute to the radiative transfer, and their weights and coupling depend on the ratio of the time step to the local particle’s mean collision time. The UGKWP method is also extended to the coupled radiation-material system. With the inclusion of energy exchange, the UGKWP method can give excellent simulation results in different regimes. A few benchmark problems are tested to show the performance of the current scheme. The accuracy and efficiency of the UGKWP method are fully confirmed. For the multiscale transport, in the cases with the co-existence of different regimes the UGKWP may improve the efficiency on several-order-of-magnitude in comparison with the purely particle methods. The UGKWP takes also advantages of both particle and macroscopic solver. The UGKWP has totally removed the ray effect and and has a much improved efficiency in comparison with DOM-type UGKS. Based on the direct modeling of the transport physics in the time step scale [30], such as the relationship between $t_{c}$ and $\Delta t$ in the current study, we have three versions of UGKS for the multiscale transport simulations, i.e., the DOM or DVM type, the purely particle formulation, and the wave-particle one. The advantages and disadvantages of these different discretization under the same UGKS framework will be further investigated, and the extension of these schemes to complex systems, such as radiative-hydrodynamics, plasma, multi-component system, and reactive flow, will be constructed.

Acknowledgement

The current research is supported by Hong Kong research grant council (16206617), National Science Foundation of China (11472219,11772281,91852114,11771035, U1530401).