Unified Gaskinetic WaveParticle Methods III: Multiscale Photon Transport
Abstract
In this paper, we extend the unified gaskinetic waveparticle (UGKWP) method to the multiscale photon transport. In this method, the photon free streaming and scattering processes are treated in an unsplitting 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 gaskinetic 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 wavedescribed and particledescribed photons are automatically adapted according to the numerical resolution and transport regime. Compared to the 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 type UGKS and IMC method, the UGKWP method may have severalorderofmagnitude reduction in computational cost and memory requirement in solving some multsicale transport problems.
keywords:
Radiative transfer equations, Diffusion equation, Waveparticle formulation, Monte Carlo particle method, Unified gaskinetic scheme1 Introduction
The radiative transfer equation describes photon propagation in the background medium and has important applications in the fields of astrophysics davis2012 (), atmospheric physics marshak20053d (), optical imaging klose2002optical () and so on. In this paper, we focus on the gray radiative transfer equation with isotropic scattering, which reads
(1) 
where is the specific intensity in time , space , and angle . Here is the speed of light, is the scattering coefficient, is the absorption coefficient, and 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) hunter2013comparison (); coelho2014advances (); chen2015 (); roos2016conservation () and the moment methods frank2006partial (); carrillo2008numerical (); vikas2013radiation (); alldredge2016approximating (). The second category is the stochastic methods, for example, the Monte Carlo method fleck1971 (); lucy1999computing (); hayakawa2007coupled (). 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 fleck1984random (); giorla1987random (); densmore2007hybrid (); densmore2012hybrid (), 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 asymptoticpreserving (AP) DOM klar1998asymptotic (); naldi1998numerical (); jin1999efficient (); jin2000uniformly (); mieussens2013asymptotic (); sun2015asymptotic1 (); sun2015asymptotic (); sun2017multidimensional (); sun2017implicit (); sun2018asymptotic () . One of the examples is the unified gaskinetic 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 xu2010unified (); xu2015direct (); jiang () and has been applied to the field of radiative transfer mieussens2013asymptotic (); sun2015asymptotic1 (); sun2015asymptotic (); sun2017multidimensional (); sun2017implicit (); sun2018asymptotic (), plasma transport liu2017unified (), and disperse multiphase flow LIU2019264 (). 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 xu2015direct (); xuliupof (). 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 timedependent 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 xu2010unified (). However, the corresponding purely particle version and waveparticle version of UGKS can be constructed as well liuzhuxu (); zhuliuzhongxu (). In this work, we develop a novel unified gaskinetic waveparticle (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 particlebased unified gaskinetic 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 freestream 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 resampled 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 xu2015direct () 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 gaskinetic 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 radiationmaterial 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 gaskinetic scheme (UGKS) was initially developed for the problems in the field of rarefied gas dynamics xu2010unified (); xu2015direct (), and has also been successfully applied to problems in radiative transfer under the finite volume framework mieussens2013asymptotic (); sun2015asymptotic1 (); sun2015asymptotic (); sun2017multidimensional (); sun2017implicit (); sun2018asymptotic (). 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
(2) 
which gives a nondimensional equation
(3) 
where . The nondimensionalization is used same as that in mieussens2013asymptotic (), and 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 . Denote as the area covered by this cell, as the interfaces of its boundary , as the outward normal of , and as the volume of cell . Define
to be the averaged specific intensity over the spatial cell, and
to be the averaged energy density function 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,
(4) 
and
(5) 
where the microscopic and macroscopic flux terms are respectively
(6) 
and
(7) 
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
(8) 
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 . 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 and at the beginning of each time step. The discretization of is based on the DOM in the original UGKS. It has been proved in mieussens2013asymptotic () that when equals , 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
In the diffusion limit, with a uniform mesh the UGKS scheme becomes a standard central difference method for the limit diffusion equation as tends to . Specifically, for 1D case with uniform grid, in the diffusive limit the scheme is
More details on the asymptotic analysis of the UGKS for the radiative transfer equation can be found in mieussens2013asymptotic () and sun2015asymptotic1 ().
Following the methodology of the UGKS, we will construct two particlebased algorithms with multiscale transport property for recovering transport physics from the kinetic scale to the hydrodynamic scale, i.e., the unified gaskinetic particle (UGKP) method and the unified gaskinetic waveparticle (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 resample these particles subsequently. The multiscale particle methods for equation (8) are constructed through the tracking and resampling of particles with the help of updated macroscopic variables. We will first introduce the UGKP method in the next section.
3 The unified gaskinetic particle (UGKP) method
In this section, we present a unified gaskinetic 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,
(9) 
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 twodimensional case on uniform Cartesian mesh. Its extension to nonuniform 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),
(10) 
The local equilibrium can be expanded as
(11) 
and the integral solution can be written as
(12) 
The first order approximation gives
(13) 
and the second order approximation gives
(14) 
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 to free stream and has a probability of to get scattered. The scattered particle follows a velocity distribution , 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 resampled at from their evolved distribution, as shown in Fig. 2. The freestreaming particles are named as the particles, and the scattered particles are named as the particles. The total energy of the particles is called 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 is defined as the time when the first scattering happens to the particle. According to Eq.(12), the time can be obtained for each particle as
where is a random number which follows a uniform distribution on . Simulation particles can be categorized into two groups based on . The first group of simulation particle satisfies
and does not suffer any collision during a time step , and this particle position is updated by
The second group of particle satisfies , those particles suffer one or multiple times of scattering during , with the assumption for simplification. Therefore, these particles will only be accurately tracked for and be resampled from at . 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
(15) 
where the flux is the diffusive flux and is the contribution of the free streaming flux. Based on (12), the diffusive flux is computed same as that in the UGKS,
We adopt the piecewise linear reconstruction for ,
(16) 
and the explicit central difference discretization of ,
where , are the barycenters of cells and with the same interface . Then the diffusive flux is given by
(17) 
where
The free transport flux is contributed by the particle free streaming, which is computed by
(18) 
where and are respectively the weight and velocity of the particle , and is the set of all particles which are transported across interface during the free streaming process. It recovers the other part of UGKS flux
(19) 
With the diffusive flux Eq. (17) and the free transport flux Eq. (18), the macroscopic energy can be updated by
(20) 
With the updated , the energy of the scattered particles can be obtained and be used to regenerate particles,
(21) 
where is the set of particles in cell at . In summary, the evolution of particle provides the freestream flux for macroscopic energy evolution, and the updated energy 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 gaskinetic waveparticle (UGKWP) method
4.1 The multiscale evolution of particles and macroscopic variable
In UGKP, the particles are resampled from an equilibrium distribution at . According to the integral solution (10), particles have a probability to get scattered in the next time step from to . We define the scattered particles in as particles and the nonscattering particles as particles. As the particles will be deleted after within , and their contribution to the free transport flux could be computed analytically, only the particles need to be resampled and they keep on free transport in the step from to . For the particles, we only need to store the distribution function instead of resampling particles from it. Based on this observation, we propose a more efficient unified gaskinetic waveparticle (UGKWP) method.
The unified gaskinetic waveparticle(UGKWP) method improves UGKP in two aspects:

Only particles are sampled as shown in Fig. 5;

The free transport flux contributed by can be calculated analytically.
The finite volume energy evolution follows
(22) 
where the diffusive flux is computed using Eq. (17), and the free transport flux contributed by is calculated by
The free transport flux contributed from the particles in Eq. (22) is computed by particle tracking,
(23) 
where is the set of all particles which are transported across the cell interface during the free streaming process. In the particle free transport process, the resampled particles will always have and keep free transport without collision in the whole time step. The other particles from will transport according to their individual generated . With the solution of , Eq.(21) is used to get , which is divided into
and
and only particles from 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 gaskinetic waveparticle method satisfies the following properties:

The microscopic and macroscopic evolutions are consistently evaluated based on the waveparticle decomposition. The macroscopic energy density is the sum of the energy from the unscattering particles and the scattering one . The evolution of particles can be evaluated analytically.

In the diffusive limit, and , essentially no particle can be freely transported and survived within a time step. Therefore, no particle is resampled from the energy. The algorithm becomes a timeexplicit central difference solver of the diffusion equation. For instance, in 1D case with uniform grid, the scheme tends to the following limiting equation
Note that in the diffusion regime, different from DOM method there is no discrete velocity space discretization in UGKWP.

In the free transport limit, and , 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,
(24) 
Define and , then the second equation can be rewritten as
(25) 
The implicit Monte Carlo method proposed by Fleck and Cummings in fleck1971 () has been shown to be an effective technique for solving nonlinear, timedependent, 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 semiimplicit discretization for material temperature evolution will be employed. Specifically, Eq. (25) is discretized by
which gives
(26) 
With the definition
substituting Eq. (26) into Eq. (24) yields
An operator splitting approach is used to solve the above system subsequently by the linear kinetic equation
(27) 
the radiation energy exchange,
(28) 
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
(29) 
where 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 energy. Afterwards, the energy change of particles is summed up, and the material temperature and are updated by energy conservation. The above UGKWP can be further improved if the technique in sun2015asymptotic () 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 .
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 i78700K CPU and GB memory. In all simulations, the time step is determined by , with for 1D examples and 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 nondimensional linear equation
defined on the semiinfinite spatial domain with an isotropic inflow condition imposed on the left boundary. For the numerical simulation, the spatial domain is taken to be . Inflow boundary condition is imposed at with the incoming specific intensity . The initial value is for all .
The results of both the UGKWP method and the Monte Carlo method are obtained using 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 (), the intermediate (), and the diffusive () regimes, as defined in jin1998diffusive (). When 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 resampling, fewer particles are used in UGKWP than the Monte Carlo method when 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 , while Table 2 compares the maximum number of particles used. Notice that no particle is generated in UGKWP method when and the method only solves a diffusion equation. For , the Monte Carlo method takes more than half an hour while the UGKWP method takes only seconds, showing that the UGKWP method is much more efficient than the Monte Carlo method in the diffusive regimes.
MC  UGKWP  

s  s  
s  s  
s  s 
MC  UGKWP  

In Figure 7 the numerical results of the UGKWP method are compared with the solutions of the Monte Carlo method for different .
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 , 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 setup 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
Both the UGKWP method and the Monte Carlo method use grids in space.
We calculate two test cases. For the first test case, we take . 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 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 particles and takes seconds. The UGKWP method uses a maxium of particles and takes 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 , 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 and they are almost identical except for the statistical noise. The Monte Carlo method uses a maximum of particles and spends seconds while the UGKWP method uses a maximum of particles and spends 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
for semiinfinite domain with a computational one . We take absorption coefficient to be , the speed of light to be , the parameter and the specific heat to be . The initial material temperature is set to be and initially material and radiation energy are at equilibrium. A constant isotropic inflow radiation intensity with an equivalent radiation temperature of is imposed on the left boundary. Both the UGKWP method and the Monte Carlo method use grids in space. For this example, the implicit Monte Carlo method uses a maximum of particles, while the UGKWP method uses a maximum of particles. It takes the implicit Monte Carlo method seconds to compute until and the UGKWP method seconds.
In Figure 9(a), the computed radiation wave front at time , and 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 Linesource problem in purely scattering homogeneous medium
The next problem we look at is the linesource problem in purely scattering medium. In this test case, we consider equation (3) for and . The spatial 2D computational domain is . The initial density distribution is
which means that initially all particles are concentrated at and they spread out over the time. This is a particularly difficult problem for the based methods because they suffer from severe ray effect. Previous studies have used this problem for comparisons between methods brunner2002forms () and as a test case for schemes to mitigate the ray effect camminady2019ray ().
Both the UGKWP method and the Monte Carlo method use cells in space. The UGKWP method uses particles while the Monte Carlo method uses 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 seconds for the UGKWP method and seconds for the Monte Carlo method. The UGKWP method is slightly more expensive due to its additional overheads.
Fig.10 compares the solution of between the UGKWP method, the Monte Carlo method and the method for time . For the computation, we use an equalweight quadrature set, cells for the spatial mesh and the DOMUGKS. 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 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 mcclarren2010robust () by considering the following equation
The diffusion term is solved by using a backward Euler method. The twodimensional physical space 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 and ; the red regions are pure absorbers with and . In the white region in the center there is an isotropic source . The initial material temperature is and initially the radiation and material temperature are at equilibrium. We take , and to simulate the case.
Both the UGKWP method and the Monte Carlo method use a mesh of in physical space. The results at time 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 particles and took seconds, while the UGKWP method uses particles and took 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 mcclarren2010robust (). We study equation (24) with . 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 . Initially radiation and material temperature are at equilibrium and the initial material temperature is .
The UGKWP method uses grids in space while the Monte Carlo method use grid in space. We compute for . The Monte Carlo method uses a maximum of particles and takes hours. The UGKWP method uses a maximum of particles and takes seconds, therefore the UGKWP method is 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 mcclarren2010robust () pointed out solution to this problem should contain two characteristics: the central block should be heated nonuniformly; 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 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 gaskinetic waveparticle 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 radiationmaterial 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 coexistence of different regimes the UGKWP may improve the efficiency on severalorderofmagnitude 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 DOMtype UGKS. Based on the direct modeling of the transport physics in the time step scale xu2015direct (), such as the relationship between and 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 waveparticle 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 radiativehydrodynamics, plasma, multicomponent 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).
References
 (1) S. W. Davis, J. M. Stone, Y.F. Jiang, A radiation transfer solver for athena using short characteristics, The Astrophysical Journal Supplement Series 199 (1) (2012) 9.
 (2) A. Marshak, A. Davis, 3D radiative transfer in cloudy atmospheres, Springer Science & Business Media, 2005.
 (3) A. D. Klose, U. Netz, J. Beuthan, A. H. Hielscher, Optical tomography using the timeindependent equation of radiative transferâpart 1: forward model, Journal of Quantitative Spectroscopy and Radiative Transfer 72 (5) (2002) 691–713.
 (4) B. Hunter, Z. Guo, Comparison of quadrature schemes in dom for anisotropic scattering radiative transfer analysis, Numerical Heat Transfer, Part B: Fundamentals 63 (6) (2013) 485–507.
 (5) P. J. Coelho, Advances in the discrete ordinates and finite volume methods for the solution of radiative heat transfer problems in participating media, Journal of Quantitative Spectroscopy and Radiative Transfer 145 (2014) 121–146.
 (6) S.S. Chen, B.W. Li, Y.S. Sun, Chebyshev collocation spectral method for solving radiative transfer with the modified discrete ordinates formulations, International Journal of Heat and Mass Transfer 88 (2015) 388–397.
 (7) T. Roos, T. Harms, C. du Toit, Conservation of scattered energy and asymmetry factor in the new rotationally symmetric spherical discretisation scheme, International Journal of Heat and Mass Transfer 101 (2016) 205–225.
 (8) M. Frank, B. Dubroca, A. Klar, Partial moment entropy approximation to radiative heat transfer, Journal of Computational Physics 218 (1) (2006) 1–18.
 (9) J. A. Carrillo, T. Goudon, P. Lafitte, F. Vecil, Numerical schemes of diffusion asymptotics and moment closures for kinetic equations, Journal of Scientific Computing 36 (1) (2008) 113–149.
 (10) V. Vikas, C. Hauck, Z. Wang, R. O. Fox, Radiation transport modeling using extended quadrature method of moments, Journal of Computational Physics 246 (2013) 221–241.
 (11) G. W. Alldredge, R. Li, W. Li, Approximating the method by the extended quadrature method of moments for radiative transfer in slab geometry, Kinetic & Related Models 9 (2).
 (12) J. Fleck, J. Cummings, An implicit Monte Carlo scheme for calculating time and frequency dependent nonlinear radiation transport, J. Comput. Phys. 8 (3) (1971) 313–342.
 (13) L. Lucy, Computing radiative equilibria with Monte Carlo techniques, Astronomy and Astrophysics 344 (1999) 282–288.
 (14) C. K. Hayakawa, J. Spanier, V. Venugopalan, Coupled forwardadjoint Monte Carlo simulations of radiative transport for the study of optical probe design in heterogeneous tissues, SIAM Journal on Applied Mathematics 68 (1) (2007) 253–270.
 (15) J. Fleck Jr, E. Canfield, A random walk procedure for improving the computational efficiency of the implicit Monte Carlo method for nonlinear radiation transport, Journal of Computational Physics 54 (3) (1984) 508–523.
 (16) J. Giorla, R. Sentis, A random walk method for solving radiative transfer equations, Journal of Computational Physics 70 (1) (1987) 145–165.
 (17) J. D. Densmore, T. J. Urbatsch, T. M. Evans, M. W. Buksas, A hybrid transportdiffusion method for Monte Carlo radiativetransfer simulations, Journal of Computational Physics 222 (2) (2007) 485–503.
 (18) J. D. Densmore, K. G. Thompson, T. J. Urbatsch, A hybrid transportdiffusion Monte Carlo method for frequencydependent radiativetransfer simulations, Journal of Computational Physics 231 (20) (2012) 6924–6934.
 (19) A. Klar, An asymptoticinduced scheme for nonstationary transport equations in the diffusive limit, SIAM journal on numerical analysis 35 (3) (1998) 1073–1094.
 (20) G. Naldi, L. Pareschi, Numerical schemes for kinetic equations in diffusive regimes, Applied mathematics letters 11 (2) (1998) 29–35.
 (21) S. Jin, Efficient asymptoticpreserving (AP) schemes for some multiscale kinetic equations, SIAM Journal on Scientific Computing 21 (2) (1999) 441–454.
 (22) S. Jin, L. Pareschi, G. Toscani, Uniformly accurate diffusive relaxation schemes for multiscale transport equations, SIAM Journal on Numerical Analysis 38 (3) (2000) 913–936.
 (23) L. Mieussens, On the asymptotic preserving property of the unified gas kinetic scheme for the diffusion limit of linear kinetic models, Journal of Computational Physics 253 (2013) 138–156.
 (24) W. Sun, S. Jiang, K. Xu, An asymptotic preserving unified gas kinetic scheme for gray radiative transfer equations, Journal of Computational Physics 285 (2015) 265–279.
 (25) W. Sun, S. Jiang, K. Xu, S. Li, An asymptotic preserving unified gas kinetic scheme for frequencydependent radiative transfer equations, Journal of Computational Physics 302 (2015) 222–238.
 (26) W. Sun, S. Jiang, K. Xu, A multidimensional unified gaskinetic scheme for radiative transfer equations on unstructured mesh, Journal of Computational Physics 351 (2017) 455–472.
 (27) W. Sun, S. Jiang, K. Xu, An implicit unified gas kinetic scheme for radiative transfer with equilibrium and nonequilibrium diffusive limits, Communications in Computational Physics 22 (4) (2017) 889–912.
 (28) W. Sun, S. Jiang, K. Xu, An asymptotic preserving implicit unified gas kinetic scheme for frequencydependent radiative transfer equations, International Journal of Numerical Analysis & Modeling 15.
 (29) K. Xu, J.C. Huang, A unified gaskinetic scheme for continuum and rarefied flows, Journal of Computational Physics 229 (20) (2010) 7747–7764.
 (30) K. Xu, Direct modeling for computational fluid dynamics: construction and application of unified gaskinetic schemes, World Scientific, 2015.
 (31) D. Jiang, M. Mao, J. Li, X. Deng, An implicit parallel UGKS solver for flows covering various regimes, Advances in Aerodynamics, (2019) 1:8 https://doi.org/10.1186/s4277401900085.
 (32) C. Liu, K. Xu, A unified gas kinetic scheme for continuum and rarefied flows v: multiscale and multicomponent plasma transport, Communications in Computational Physics 22 (5) (2017) 1175–1223.
 (33) C. Liu, Z. Wang, K. Xu, A unified gaskinetic scheme for continuum and rarefied flows vi: Dilute disperse gasparticle multiphase system, Journal of Computational Physics 386 (2019) 264 – 295.
 (34) K. Xu, C. Liu, A paradigm for modeling and computation of gas dynamics, Physics of Fluids 29 (026101).
 (35) C. Liu, Y. Zhu, K. Xu, Unified gaskinetic waveparticle methods i: continuum and rarefied gas flow, arXiv preprint arXiv:1811.07141.
 (36) Y. Zhu, C. Liu, C. Zhong, K. Xu, Unified gaskinetic waveparticle methods ii: multiscale simulation on unstructured mesh, arXiv preprint arXiv:1903.11861.
 (37) S. Jin, L. Pareschi, G. Toscani, Diffusive relaxation schemes for multiscale discretevelocity kinetic equations, SIAM Journal on Numerical Analysis 35 (6) (1998) 2405–2439.
 (38) T. A. Brunner, Forms of approximate radiation transport, SAND20021778, Sandia National Laboratory.
 (39) T. Camminady, M. Frank, K. Küpper, J. Kusch, Ray effect mitigation for the discrete ordinates method through quadrature rotation, Journal of Computational Physics.
 (40) R. G. McClarren, C. D. Hauck, Robust and accurate filtered spherical harmonics expansions for radiative transfer, Journal of Computational Physics 229 (16) (2010) 5597–5614.