Radiation Transfer Modeling in Atmospheric Aerosol Media Using a GPU-Implemented Lattice Boltzmann Model

1. Introduction

The aerosols are abundant in the troposphere, especially in the atmospheric boundary layer. They are either generated from natural processes (dust, sea salt, fog, cloud) or from human activities (smoke, smog, and other air pollutants). Since the aerosols are effective radiative scatter, they play an important role in atmospheric radiative energy transfer. The aerosol particle sizes are much larger than the air molecules, usually in the range of a few wavelengths of visible light to several micrometers, and the light interaction with aerosol does not follow the isotropic Rayleigh scattering law of molecular scattering. The scattering of light by aerosols is highly anisotropic and concentrated in forward or backward directions, depending on the properties and morphology of aerosols. Accurate modeling of radiation transfer in the aerosol medium is important in analysis of natural radiative energy transfer processes in atmosphere, light signal propagation, directed energy transfer, and remote sensing of atmospheric environment.

The radiative transfer process can be modeled using radiative transfer equation (RTE). Since analytical solutions are possible only for simplest, idealized conditions, different computation methods, such as discrete ordinate method (DOM) [1,2,3], finite volume method [4,5], and Monte Carlo method (MCM) [6] have been developed to solve the RTE numerically. All those numerical methods of solving RTE demand a large amount of computation time. MCM is often used as a benchmark method to validate the numerical solution of other methods. Although MCM is accurate and versatile in solving the RTE, it demands most intensive computing compared with other methods. As computation domain is getting larger or resolution of grid is finer, computation intensity grows rapidly.

An accurate and computationally efficient method is desired for modeling radiative energy transfer in the atmosphere. During recent years, great progress has been made in developing lattice Boltzmann method (LBM) for rapid computation of radiative energy transfer problems. The LBM was first discovered and developed in the fluid mechanics community [7,8,9,10,11,12,13,14,15,16] and has become one of the most effective methods in terms accuracy and computation speed for fluid flow and heat transfer simulations. LBM is based on the kinetic theory of statistical mechanics and solves the Boltzmann equation that governs the probability distribution of fluid particles. The LBM solves the Boltzmann equation for a particle at each grid point by performing collision and propagation calculations of the particle’s probability distribution function (PDF) over a discrete and symmetric lattice mesh with certain fixed directions. The macroscopic variables of fluid flow, such as u,v,w, can be derived from the statistical moments of PDF solved by LBM. The LBM flow models, implemented on a GPU computation platform, can drastically increase computation speed by about 300 times [16] compared with a CPU implementation. The LBM algorithm is well suited on GPU because LBM calculations are highly localized, data parallel. A particle collision and streaming steps in LBM only use information from its immediate neighbor grids. The RTE equation comprises terms of energy advection (propagation) and diffusion (scattering/absorption), it is natural to explore the LBM as a computational method to accelerate the computation speed in modeling radiative transfer. The LBM can be considered a direct discretization of the Boltzmann equation [7]. Because of the similarity between the RTE and Boltzmann equation [17], earlier research and development of the LBM for the RTE started with direct discretization of the RTE with respect to space, time, and angular direction. The earliest work of solving the RTE using the LBM by Giest et al. [18] was for rendering computer graphics considering the radiation effects. Asinari et al. and Mishra et al. [19,20] developed a two-dimensional (2-D) LBM for radiative transfer modeling in a particular medium. Ma et al. [21] derived an LBM used for a one-dimensional (1-D) radiation problem that compared well with an analytical solution. Bindra and Patil [22] and McCulloch and Bindra [23] also developed a 2-D LBM for the RTE for a simulation of the conjugated radiative and convective heat transfer problem. A general review on modeling neutron and photon transport using LBM is provided in [24]. The LBM was also used in a non-equilibrium radiation transfer problem [25]. Zhang et al. [26] and Yi et al. [27] derived a 2-D LBM using the Chapman–Enskog expansion for a steady-state radiative transfer problem that can deal with both thin and high optical depths. For better treatment of the radiation source term, a multi-relaxation time LBM was developed by Liu et al. [28]. Liu et al. [29] also developed an unstructured grid method in LBM for speed up computation and easier treatment boundary in a homogeneous medium. McHardy et al. [30,31] developed a 3-D LBM model using a direct discretization of the RTE and the model produced accurate results for the ballistic radiation condition in which the medium scattering albedo is less than 0.7. Mink et al. [32,33] developed a 3-D LBM method for high optical thickness situations based on the Chapman–Enskog expansion and a steady-state RTE was approximated by the Helmholtz equation and solved with the LBM. Wang et al. [34] developed a GPU accelerated LBM model for radiation transfer in clean atmosphere with molecular scattering and the great speed acceleration was realized.

Radiative transfer modeling in different aerosol mediums is more complex and less developed compared with its counterparts of radiative transfer in atmospheric molecules. One aspect of difficulty is that radiative parameters for different species of aerosols are not well understood and much more related to material properties and morphology of the aerosols. Another aspect of difficulty is the computation of the radiation transfer process due to the asymmetry scattering parameters. An LBM was developed by McHardy et al. [30] for case of Mie scattering and the results compared well with the MCM. The key signature in solving RTE using LBM for anisotropic scattering aerosol medium is that an anisotropic scattering function (SPF) must be applied in each computation grid. For anisotropic aerosol scattering of larger aerosol particles, the SPF is not uniform and must be specified with incoming angles. Furthermore, the SPE must be normalized to obey the energy conservation in the system [3,30,35,36] . SPE of aerosols can be derived from theoretical analysis of solutions of the Maxwell equations [37,38,39] or laboratory observations [40,41].

The objective of this research is to develop and validate a radiation transfer lattice Boltzmann model (RTLBM) for atmospheric aerosol medium to speed up the computation. This paper is organized as follows: section 2 describes the basic theory of the radiation transfer and derivation of the RTLBM. Section 3 shows the validation of the RTLBM results with a radiation transfer Monte Carlo model (RTMCM) to test the accuracy of the RTLBM, the results of computation times are compared for different sizes of domains to demonstrate the computation efficiency of which RTLBM is implemented on a GPU platform. The final section summaries this research and conclusion.

2. Methods

2.1. The Lattice Boltzmann Equation for the RTE of Aerosol Medium

The RTE is an energy conservation equation for radiation transfer and can be written in the following form:

$${\displaystyle \frac{1}{c}}{\displaystyle \frac{\partial R}{\partial t}}+\mathit{n}·\nabla R=-\left({\mu }_a+{\mu }_s\right)R+{\mu }_a{R}_b+{\displaystyle \frac{{\mu }_s}{4\pi }}{\int }_{4\pi }\Phi \left({\mathit{n}}^{\mathit{\text{'}}},\mathit{n}\right)R^{\text{'}}d{\Omega }^{\text{'}}$$

(1)

where R(x,$\mathit{n}$,t) is the spectral radiance at spatial point x, $\mathit{n}$ is a direction vector, c is the speed of light, ${\mu }_a\mathrm{a}\mathrm{n}\mathrm{d}{\mu }_s$ are the medium absorption and scattering coefficients, respectively. $R_b$ is the spectral radiance of blackbody emission, $\Omega$ is the solid angle. Since the scattering from aerosol particles can be treated as monochromatic, the spectral dependency of radiance can be dropped out. The $\Phi ({\mathit{n}}^{\mathit{\text{'}}},\mathit{n})$ is the SPF to direction $\mathit{n}$ from direction ${\mathit{n}}^{\mathit{\text{'}}}$. In the case of isotropic molecular scattering medium, $\Phi \left({\mathit{n}}^{\mathit{\text{'}}},\mathit{n}\right)=1.$ But for aerosol medium, the scattering is anisotropic, the scattering phase function is different from uniform, it needs an asymmetry parameter, $g$, to account for the distribution of scattering with respect to scattering angle (see Section 2.2).

In the derivation of the Boltzmann type of radiation transfer equation, the relation between the PDF at a direction i ($f_i\left(\mathit{x},t\right)$) of a virtual particle crossing a solid angle cross-section and the radiance is expressed as [33]

$$f_i\left(\mathit{x},t\right)=w_{i}R(\mathit{x},n_i,t),$$

(2)

where $w_i$ are the solid angle weights that correspond to the lattice directions. In this article, lattice model with 26 directions is used and their values are

$$w_i=\left\{\begin{array}{c}{\displaystyle \frac{1}{21}},i=1,\dots ,6\\ {\displaystyle \frac{4}{105}},i=7,\dots ,18\\ {\displaystyle \frac{9}{280}},i=19,..,26\end{array}\right.,$$

(3)

where the i directions in lattice structure are shown in Figure 1. By neglecting the blackbody emission of the aerosol particles, and substituting the absorption and scattering coefficient with extinction coefficient (or optical depth), $\beta$, and scattering albedo, α, the RTE is transformed into following Boltzmann type equation:

$${\displaystyle \frac{\partial f_i}{\partial t}}+n_i.\nabla f_i=-c\beta \left[f_i-{\displaystyle \frac{\mathsf{\alpha }}{4\pi }}{\displaystyle {\int }_{4\pi }}\Phi (n_j^{\text{'}},n_i)f_i^{\text{'}}d{\mathsf{\Omega }}^{\text{'}}\right],$$

(4)

where $\beta =$ ${\mu }_a+{\mu }_s$ , $\alpha ={\mu }_s/\beta$.

To keep the model non-dimensional for the comparisons and applications, the medium’s scattering albedo, α, and optical depth, $\beta$ , both are non-dimensional parameters, are used instead of the coefficients of absorption and scattering. The optical depth can be expressed in a normalized form

$$\beta =\left({\mu }_a+{\mu }_s\right)l_{phy},$$

(5)

where $l_{phy}=1$is a modeled normalized physical domain length. The main complexity of solving the RTE is in the scattering integral term, especially in the anisotropic scattering medium. The scattering integral is replaced with discrete form summed over all the lattice directions

$${\displaystyle \frac{1}{4\pi }}{\displaystyle {\int }_{4\pi }}\Phi \left(n_j^{\text{'}},n_i\right)f_i^{\text{'}}d{\mathsf{\Omega }}^{\text{'}}=w_i{\displaystyle \sum _{j=1}^{26}}{\Phi }_{ij}{f}_j=f_i^{eq},$$

(6)

where $f_i^{eq}$ is an equilibrium distribution of photon kinetics of particle collision. ${\Phi }_{ij}$ is a scattering matrix which represents the probability that a photon is scattered from the i to j direction, $w_i$are the solid angle weights for the integration. For an isotropic scatter particle such molecules in atmosphere, ${\Phi }_{ij}=1.$ $f_i^{eq}$ has a different mechanism in radiative transfer than in fluid flow. The equilibrium is the Maxwell and Boltzmann distribution in fluid modeling, representing the approaching fluid particle equilibrium PDF after the collision; The$f_i^{eq}$ in radiative transfer, however, represents the equilibrium PDF of photon after scattering by the medium particle.

The Equation (4) can be discretized in time, space, and velocity, the resulting lattice Boltzmann equation is as follows [7]:

$$f_i(\mathit{x}+c_i\mathrm{\Delta }t,t+\mathrm{\Delta }t)-f_i(\mathit{x},t)=-c_i\mathrm{\Delta }t\beta f_i+c_i\mathrm{\Delta }t\alpha \beta f_i^{eq},$$

(7)

where $c_i$is the discrete lattice velocity in different directions of a lattice:

$$c_i=\left\{\begin{array}{c}\left(\pm \mathrm{1,0},0\right)\left(0,\pm \mathrm{1,0}\right)\left(\mathrm{0,0},\pm 1\right)i=1,\dots ,6\\ \left(\pm 1,\pm 1,0\right)\left(\pm \mathrm{1,0},\pm 1\right)\left(0,\pm 1,\pm 1\right)i=7,\dots ,18\\ \left(\pm 1,\pm 1,\pm 1\right)i=19,\dots ,26\end{array}\right.$$

(8)

An alternative explanation and derivation of LBM for radiative transfer problem start with the microscopic Boltzmann transport equation. By using the Chapman-Enskog expansion, the macroscopic RTE equation can be recovered [26,27,28,29]. The Chapman-Enskog expansion analysis is fundamental theory that connects the microscopic particle kinetics to the collective macroscopic transport equations. Detailed description about theory and applications of the LBM in general is referred to monograph [13].

After many time steps of collision and streaming iterations, the solution of the equation approaches a steady solution of error with an error tolerance. The radiation intensity, I(x,t), can be computed from the first order statistical moment of the photon PDF. An intensity of radiation field, $I\left(\mathit{x},t\right),$ is defined as

$$I\left(\mathit{x},t\right)={\displaystyle {\sum }_{i=1}^{26}}{f}_i\left(\mathit{x},t\right)$$

(9)

2.2. Normalization of Aerosol Scattering Phase Functions

The scattering of light waves by spheric particles that are comparable in size to the wavelength of the light is described by Mie theory. The scattered light is more intense in forward or backward directions, and the intensity is dependent on wavelength of the light, property, and size of the particle. To avoid much of the complexity of Mie scattering phase function, an approximate phase function $\Phi$ of Henyey-Greenstein (HG) is used. The HG function is expressed as

$$\Phi \left(\theta \right)={\displaystyle \frac{1-g^2}{{\left[1+g^2-2g\mathrm{c}\mathrm{o}\mathrm{s}\left(\theta \right)\right]}^{3/2}}},$$

(10)

where $\theta$ is the angle from incoming radiation to scattering direction, $g$ is the asymmetry factor. To conserve the scattered energy, the HG scattering phase function needs to be normalized to assure the scattering energy to be conserved with the incoming energy and overall asymmetry factor is unaltered [35,36]:

$${\displaystyle \frac{1}{4\pi }}{\displaystyle {\sum }_{j=1}^N}{\Phi }_{ij}{w}_j=1i=1,\dots ,26,$$

(11)

$${\displaystyle \frac{1}{4\pi }}{\displaystyle {\sum }_{j=1}^N}{\Phi }_{ij}{w}_j\cos {\theta }_{ij}=gi=1,\dots ,26,$$

(12)

where the subscripts i and j represent the discrete directions i and j, N=26 of the discrete lattice directions, $w_j$ are the solid angle weights that correspond to the lattice directions. McHardy et al. [30] has used an older version of Hunter and Guo [35], which is more complex and computationally intensive for which solves an under-determined liner system involving many interactions. In this research, we chose a new method developed by Hunter and Guo [36], an explicit system which is much simpler than their previous algorithm. The new algorithm was also tested to its accuracy compared with the older algorithm [36]. The new normalization procedure the forward (f) and backward (b) scattering components was altered according to the following formulas,

$${\tilde{\Phi }}_{if}=\left(1+A_f\right){\Phi }_{if}i=1,\dots ,26,$$

(13)

$${\tilde{\Phi }}_{ib}=\left(1+B_b\right){\Phi }_{ib}i=1,\dots ,26.$$

(14)

This normalization can conserve the scattered energy while simultaneously maintaining the original value of asymmetry factor. Normalization vector parameters$A_{f}and{B}_b$ are computed as follows :

$$A_f={\displaystyle \frac{1}{2{\Phi }_{if}{w}_f}}\left[4\pi \left(1+g\right)-{\displaystyle {\sum }_{j=1}^N}{\Phi }_{ij}{w}_j(1+\cos {\theta }_{ij})\right],$$

(15)

$$B_b={\displaystyle \frac{1}{2{\Phi }_{ib}{w}_b}}\left[4\pi \left(1-g\right)-{\displaystyle {\sum }_{j=1}^N}{\Phi }_{ij}{w}_j(\cos {\theta }_{ij}-1)\right].$$

(16)

The new algorithm is based on the notion that strong forward scattering HG phase function is significantly larger than the remaining terms [39], only the large terms of forward or backward (for mainly backward scattering aerosols) are altered for scattering energy conservation. In this study, we only consider the forward scattering aerosols. Figure 2 displays the examples of HG normalization results for the different asymmetry parameters. The results show that only the forward direction is altered in the HG function and other values of other directions are on the non-normalized HG curves.

In many situations, the aerosol particles are not spheric but rather in irregular shapes such as ellipsoids, columns, clusters of smaller particles. The laboratory test [40,41] can provide SPF, $\Phi \left(\theta \right),$for single scattering of an aerosol particle. The T-Matrix method can also be applied for single scattering of irregular shaped aerosol particles . The retrieval of SPF for mixed sized aerosol particles can also be derived from the T-Matrix method [37,38,39]. It is not feasible to compute a multiple scattering of countless aerosol particles. It is also to be emphasized that SPF also needs to be normalized as demonstrated to ensure energy and $g$ parameter conservations in the multi-scattering medium.

2.3. Boundary Conditions

Three types of boundary conditions were used in this study. The first type is the beam type radiation, represents radiation was from a specified unit beam radiation with a specified elevation angle. The second type of boundary condition is diffused incoming radiation from a boundary. The third type of boundary condition is a non-emitting black boundary condition allowing the radiation to be absorbed completely when the photo hits the boundary. The fourth type of boundary condition is a periodic boundary condition along the lateral boundary planes, representing the simulation domain is uniform out of the simulation domain. The examples of the boundary conditions can be described as follows in the RTLBM (see Figure 1 for the lattice directions):

$$\left\{\begin{array}{c}\begin{array}{c}f6=I_b,andf\mathrm{13,14,17,18,19,22,24,25}=0,forperpendicularbeamfromtop;\\ f1,\mathrm{7,9},\mathrm{11,13,19,21,23},25=I_b{\displaystyle \frac{w_i}{\overline{w}}},fordiffuedradiationfromwestboundary;\\ f2=f1,f8=f9,f10=f7,f12=f13,f14=f11,f20=f19,f26=f25,\end{array}\\ f22=f21,f24=f23,foreast-westperiodicboundaryconcition;\\ f1,\mathrm{7,9},\mathrm{11,13,19,21,23},25=0,fornon-emittingblackwestboundary.\end{array}\right.$$

(17)

where$I_b$ is the radiative intensity at the incoming boundary, and $\overline{w}$ is the sum of all weights in the incoming direction.

The RTMCM used in this study is the same as the one that was used in our previous study. For brevity, interested reader is referred to Wang et al. [34].

3. Results

To evaluate the accuracy of the newly developed LBM model for the radiative transfer in the aerosol medium, MCM results were used for this purpose. MCM [34] was also used for evaluation purposes in previous LBMRT modeling of radiative transfer in the isotropic molecular Rayleigh scattering in the atmosphere. Despite MCM huge computational power requirement, MCMs are used for benchmark validation tools in other types of model developments. We also follow the suit of this practice since the observational data is scarce. Since the radiation transfer in aerosol is mainly in high scattering diffusive regime, we focused the following simulations on high scattering albedo.

3.1. Radiation Transfer Simulations in Aerosol Medium with a Beam Incoming Boundary

In this validation of the LBM model for anisotropic radiation transfer in an aerosol medium. A cubic volume of aerosol medium with specified scattering phase function and optical depth were specified. The optical depth of the aerosol medium is 33.333, and the scattering phase function parameter is 0.85. The top boundary is incoming with beam radiation perpendicular to the top plane the computational domain with normalized unit intensity (Figure 3, left panel), and other boundaries are all assumed to be black and non-emitting.

Figure 3 (right panel) are samples of radiative intensities along the Z at the centerline (X=Y=50 grid) computed with both RTMCM and RTLBM methods. Figure 4 displays a 3D dimensional perspective visualization of the computed radiative intensity. Two vertical plane slices, one is parallel to X axis, and another is parallel to Y axis of the simulation domain are shown in the figure. The left panel is the computational results of RTLBM, and right panel is the results from the RTMCM, respectively. Since the model results are all normalized values, we use the root mean square difference (RMSD) to evaluate the average difference between the RTLBM and RTMC results. The RMSD is 0.0307 between the RTLBM and RTMCM for this case of beaming from top boundary to an aerosol medium.

3.2. Radiation Transfer Simulations in Aerosol Medium with a Diffused Incoming Boundary

In this comparison of radiation transfer simulations using the RTLBM and RTMCM methods, the aerosol medium parameters were exactly same as in case 1 stated in section 3.1. The incoming boundary condition was from west plane, and the radiation was a diffused with in RTLBM (same incoming in lattice direction of 1,7,9,11,13,19,21,23,25 (see Figure 1) . In the RTMCM, radiation was allowed to come on west plane in any direction.

McHardy et al. [30] also has done a simulation using their LBM with same setup as current study. The only difference between their implementation of RTLBM is that they used an older method of normalization in phase scattering function, which is much more complex than our normalization of phase scattering function. Figure 3 (right panel) displays the radiative intensity along the centerlines in X direction at Y=Z=0.5. Two simulations of RTLBM and another simulation with RTMCM are quite close along this sampling line.

Figure 5 displays vertical cross-sections of RTLBM and RTMCM simulation results of diffused incoming radiation process in an aerosol medium along the X axis at the Y=Z=0.5. The results from both simulations generally agreed reasonably well, but some minor differences near the incoming boundary. The RTMCM result near the incoming boundary had larger radiative intensity than that in the RTLBM simulation result. The differences between the LBM and MC models also showed similar differences, RMSD=0.0145.

Figure 5. Left panel: a sketch of the 3D computation domain, west diffused income boundary radiation and aerosol medium. Diffuse radiation was from west boundary. Medium aerosol scattering albedo=1, optical depth=33.333, and scattering phase function parameter= 0.85; Right panel: a comparison of center horizontal symmetry line samples (at Y=Z=0.5) of the RTLBM, RTMCM, and McHardy LBM [28] results.

Figure 5. Left panel: a sketch of the 3D computation domain, west diffused income boundary radiation and aerosol medium. Diffuse radiation was from west boundary. Medium aerosol scattering albedo=1, optical depth=33.333, and scattering phase function parameter= 0.85; Right panel: a comparison of center horizontal symmetry line samples (at Y=Z=0.5) of the RTLBM, RTMCM, and McHardy LBM [28] results.

Figure 6. A side by side comparison of RTLBM and RTMCM results at the center plane cuts (Y=0.5) in the radiative filed. Diffuse radiation was from west boundary. Medium aerosol scattering albeo=1, optical depth=33.333, and scattering phase function parameter= 0.85.

Figure 6. A side by side comparison of RTLBM and RTMCM results at the center plane cuts (Y=0.5) in the radiative filed. Diffuse radiation was from west boundary. Medium aerosol scattering albeo=1, optical depth=33.333, and scattering phase function parameter= 0.85.

3.3. A Case in Mixed Aerosol and Molecular Scattering Media in the Domain

One of prevalent situations for atmospheric radiation transfer is the clear atmosphere with some cloud coverage or air pollutant in some areas. This is a typical case of mixed media with Mie scattering (cloud occupied volume) and Rayleigh scattering (clean air atmosphere). In this simulation, a typical cumulus cloud area with optical depth of 20 and scattering albedo of 1, and g=0.85 was set up in the volume (x=y=0.45~0.55, z=0.2~0.3). The rest of area has a clear atmosphere with isotropic Rayleigh scattering optical depth of 0.08 and scattering albedo 0.7. The Rayleigh scattering was described in first version of the RTLBM model [34]. The identity flag was set up in the domain for tagging the computation grid occupied by the different scattering medium. The same domain was set up for the RTMCM simulation. Periodic boundary conditions were applied in lateral boundaries.

Figure 7 (left panel) shows a model domain setup as described . Figure 7 (right panel) displays a sample vertical line at center (X=Y=0.5). The RTLBM and the RTMCM produced similar results both in the Rayleigh (clean air) and Mie (cloud) scattering regions in the domain. The RTLBM simulation gave a slightly smaller value of radiative intensity along the vertical line than that of RTMCM simulation, and the location from the peak radiative intensity locations by Mie scattering had some differences. The simulation results from the RTLBM and RTMCM were also cut through the symmetric plane (Y=0.5) and displayed in Figure 8. The radiative intensity showed large values around the cloudy volume. At the cloud shadow area below the cloud showed (Z<0.2) large decrease of radiative intensity. One pattern of radiative ray effect was shown in this LBM simulation. This ray effect error was caused by the limited ray number in the LBM. The ray effect error was also evident in DOM type model due to limited number of rays. This ray effect error is more severe in the case of small and narrow medium cases [30,32]. The overall difference between the RTLBM and RTMCM was evaluated at RMSD value of 0.0210.

3.3. GPU Computation of Radiation Transfer in Aerosol Medium

One of rationale to develop the RTLBM is mainly for its computational speed and reasonably accurate results. The RTMCM model would be the best choice due to its accuracy and versatility in treatment of different types of boundary conditions. However, the RTMCM demands tremendous computational power compared with other methods because the RTMCM samples all the possible directions for billions of the photon transport processes. In all our simulations described in this section with 101³ grids, the computation time for a same single CPU takes approximately 486 times of CPU time using the RTMCM model (6 billion photons) compared with the RTLBM (Table 1). The ratio of computation time (RTMCM/RTLBM) on same CPU is even large if there is a large domain in the computation or more photons are sampled for smoother fields. In the following discussion, we will focus on comparing the computation speed of RTLBM model implemented on the CPU and GPU platforms.

As discussed in the introduction, one of the advantages of the LBM type model is the intrinsic parallelism in the modern GPU platform because the data locality is satisfied and time-stepping is explicit in the LBM algorithm. In the GPU implementation, each computational grid cell in the domain is assigned to a GPU thread. Modern GPU has several thousands of compute cores and are thus able to run thousands of threads simultaneously. The information in the grid is read, processed, and put back in the same memory by this thread. The code uses the Nvidia CUDA^® (Common Unified Architecture). The 3D computation grids are mapped to 1D memory, and the data access uses the array structure for memory coalescing. There are many thread blocks in the code, and each thread block does the same computation. Each thread block uses local fast shared memory, avoiding unnecessary local to global memory movement. In GPUs, threads execute in lockstep in sets called warps. The threads within each warp need to load memory together to use the hardware most effectively. This is called memory coalescing. In our implementation, we manage this by ensuring that threads within a warp are accessing consecutive local memory as often as possible. Current version of GPU implementation is in a single GPU processor. We are still in the process of implementing the RTLBM on a multi-GPU computation platform. Here we only report its efficiency for three cases in this paper, as indicated in Table 2. The Tesla A100 GPU (6912 cores, 80 GB memory) has more than 34 times speed-up compared to a CPU (AMD 7713 2 GHz) for a domain size of 101³ grid points. The CPU/GPU computation time ratio increases as the domain size increases, 182 times for a domain of 501*501*301, and 261 times for a domain with 501³ grid points. These results also confirmed that the LBM is a generally memory limited computation algorithm. The LBM with a GPU has shown to be more effective in numerical simulation of turbulent flow in urban environments with at least a 200 to 500 times speed-up (CPU/GPU time ratio) on early generation of the GPU type [14,15].

In the more recent GPU computation, such as AI, it is believed there are more advantages to use the easier Python implementation of the Cuda kernel for better code readability, data flow among the kernels, and rich open graphical and statistical library. The RTLBM is also implemented in Python for this advantage. As shown in Table 2 and Table 3. The speed up factors are quite similar for the C++ and Python implementation.

5. Summary and Conclusions

A fast and reasonable accurate radiative transfer model in atmospheric aerosol medium, RTLBM, was developed using the lattice Boltzmann method. A new normalization algorithm developed by Hunter and Guo [36] was used for application in this model. The RTLBM was evaluated by comparison with results from a Monte Carlo model. The test results indicated that RTLBM is reasonably accurate in the situations of direct beam and diffused incoming boundary conditions with high optical depth and scattering albedo. RTLBM also validated the situation of patch aerosol in the clean atmospheric, emulating the mixed radiation transfer media in Rayleigh and Mie scattering regimes. Like other model methods, such as DOM and FV, the RTLBM also had minor error in the situation of small of narrow patches of media due the relatively a small number of sampling directions. This error can be reduced by dense up the computational grids in the smaller region. The overall quantitate measure of difference between RTLBM with a RTMCM is that RMSD is less than 0.03 for all three test cases.

The motive for developing RTLBM is to develop RTLBM to reduce the massive computation expense compared with other numerical methods. The speedy computation of the RTLBM is based on its special algorithm of streaming and collision. As test results shown, RTLBM can easily be 260 times faster on a GPU (Nvidia A100, not even latest faster B100) than on a CPU, drastically reduced the high computation cost. This test results are comparable to the GPU implement RTLBM [34] of clear atmosphere where scattering is the isotropic Rayleigh regime. This computation speed improvement can significantly reduce the computation burden for modeling the radiative transfer process in complex atmospheric boundary layer such as urban area, where the scattering media and boundaries are complex and highly three dimensional.

Future research and development efforts will be focused on following: (1) application LBM in the plant canopy radiation modeling; (2) different boundaries and two streams (solar and terrestrial) of modeling of urban radiation energy transfer; (3) multi-GPU implementation of the RTLBM for further reduction of computation cost.