A Multigroup Moment-Accelerated Deterministic Particle Solver for Time-dependent Thermal Radiative Transfer Problems

TL;DR

This paper introduces a novel Lagrangian particle-based solver for time-dependent thermal radiative transfer that combines advantages of particle and deterministic methods, enabling accurate results with large time steps.

Contribution

It presents a new characteristic-based transport algorithm within a moment-accelerated framework, improving efficiency and robustness for TRT simulations.

Findings

Achieves accurate results with fewer particles and larger time steps.

Combines benefits of particle and grid-based methods.

Handles stiff absorption and reemission self-consistently.

Abstract

We propose an efficient, robust, Lagrangian (characteristic-based) transport solver for the time-dependent thermal radiative Transfer (TRT) applications within the context of a moment-accelerated (High-Order/Low-Order, HOLO) algorithm. This novel transport algorithm inherits the best features of both particle methods (e.g., time accuracy, phase-space adaptivity, positivity) and deterministic, grid-based methods (e.g., no stochastic noise). Particles are evolved by the method of characteristics, with a time-dependent weight that accounts for stiff absorption and reemission process self-consistently. As a result, this approach is able to obtain accurate results while employing large time steps and a moderate number of particles (compared to Implicit Monte Carlo).

Figures

Tables

Table 1. Table 1: Problem specifications for the Marshak wave problems

	Problem 1	Problem 2
	(Thick)	(Thin)
$x$ [cm]	0.25	2.0
$t_f$ [s]	$2\times {10}^{-8}$	$5\times {10}^{-8}$
$\rho$ [g/${\mathrm{cm}}^3$]	1.0	1.0
$c_v$ [erg/g-eV]	$3\times {10}^{12}$	$1.3874\times {10}^{11}$
$T_{bc}$ [eV]	1000	150
$T_0$ [eV]	0.025	0.025
$\alpha$ [${\mathrm{eV}}^3-{\mathrm{cm}}^2$/g]	$1.0\times {10}^{12}$	$1.0\times {10}^6$
$\Delta x$ [cm]	0.005	0.025

Table 2. Table 2: CPU time distribution of DP-HOLO in gray optically thick Marshak wave problem. The values in parenthesis indicates % CPU time with respect to the total CPU time.

# of particles	LO CPU time [s]	HO CPU time [s]
8	7.1 (65.6)	3.7 (34.4)
16	7.1 (48.9)	7.4 (51.1)
32	7.2 (32.6)	14.9 (67.4)
64	7.2 (19.4)	29.7 (80.6)
128	7.2 (10.8)	59.5 (89.2)
512	7.3 (2.9)	241.3 (97.1)
2048	7.4 (0.7)	858.3 (99.3)

Table 3. Table 3: CPU time distribution of DP-HOLO in gray optically thin Marshak wave problem. The values in parenthesis indicates % CPU time with respect to the total CPU time.

# of particles	LO CPU time [s]	HO CPU time [s]
8	10.5 (76.4)	3.3 (23.6)
16	10.5 (59.7)	7.1 (40.3)
32	10.7 (45.3)	12.9 (54.7)
128	10.5 (17.0)	51.4 (83.0)
512	10.7 (4.9)	206.5 (85.1)
2048	10.8 (1.3)	850.2 (98.7)

Table 4. Table 4: Problem specifications for the Marshak wave problems

$x$ [cm]	10.0
$t_f$ [s]	$1\times {10}^{-9}$
$\rho$ [g/${\mathrm{cm}}^3$]	1.0
$c_v$ [erg/g-eV]	$1\times {10}^{10}$
$T_{bc}$ [eV]	1000
$T_0$ [eV]	10.0
$\alpha$ [${\mathrm{eV}}^3-{\mathrm{cm}}^2$/g]	$1.0\times {10}^9$

Table 5. Table 5: CPU time distribution of DP-HOLO in multifrequency Larsen’s problem. The values in parenthesis indicates % CPU time with respect to the total CPU time.

# of particles	LO CPU time [s]	HO CPU time [s]
8	5.7 (54.6)	4.7 (45.4)
16	5.7 (41.5)	8.0 (58.5)
32	5.8 (28.6)	14.4 (71.4)
128	6.0 (10.0)	55.4 (90.0)
512	6.0 (2.7)	217.8 (97.3)
2048	6.1 (0.7)	867.8 (99.3)

Table 6. Table 6: Relative CPU increase for frequency-independent Marshak problem.

	$CPU\left(NG\right)/CPU\left(NG=1\right)$
NG	IMC	DP total	DP HO	DP LO
4	3.4	1.3	1.8	1.1
8	6.9	1.5	2.9	1.1
16	13.9	2.1	5.2	1.1
32	27.5	3.1	9.7	1.1
64	54.6	5.2	18.5	1.1

Table 7. Table 7: Relative CPU increase in Larsen’s problem.

	$CPU\left(NG\right)/CPU\left(NG=1\right)$
NG	IMC	DP total	DP HO	DP LO
4	61.2	1.2	1.9	1.0
8	58.4	1.8	3.5	1.3
16	89.1	2.1	5.4	1.2
32	167.5	3.2	10.0	1.4
64	338.4	5.4	19.0	1.7

Equations

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mathbf{\Omega }\cdot \nabla I+\sigma I$$

$$\rho c_v\frac{\partial T}{\partial t}+\sigma ac{T}^4-\sigma cE$$

$$B(\nu ,T)\equiv \frac{2h{\nu }^3}{c^2}\frac{1}{e^{h\nu /kT}-1}.$$

$$\frac{1}{c}\frac{\partial I}{\partial t}+\mu \frac{\partial I}{\partial x}+\sigma I$$

$$\rho c_v\frac{\partial T}{\partial t}+\sigma ac{T}^4-\sigma cE$$

$$I(x,\mu ,t)=\sum _{p^{\prime }}^{N_p}{w}_{p^{\prime }}(t)\delta (x-x_{p^{\prime }}(t))\delta (\mu -{\mu }_{p^{\prime }}(t)),$$

$$\frac{d{w}_p}{dt}+\sigma c{w}_p=cQ(x_p,t),$$

$$w_p(t)=w_p(0)e^{-{\int }_0^t\sigma cd{t}^{\prime }}+{\int }_0^t{e}^{-{\int }_{t^{\prime }}^t\sigma cd{t}^{\prime \prime }}cQ(x_p(t^{\prime }),t^{\prime })d{t}^{\prime }.$$

$$x_p(t)$$

$${\mu }_p(t)$$

$$\frac{\partial E}{\partial t}+\frac{\partial F}{\partial x}+\sigma cE=\sigma ac{T}^4,$$

$$\frac{1}{c}\frac{\partial F}{\partial t}+\frac{c}{3}\frac{\partial E}{\partial x}+\sigma F=0,$$

$$\rho c_v\frac{\partial T}{\partial t}+\sigma ac{T}^4-\sigma cE=0.$$

$$\frac{E_i^{n+1,HO}-E_i^{n,HO}}{\Delta t_n}+\frac{{\bar{F}}_{i+1/2}^{n+1,HO}-{\bar{F}}_{i-1/2}^{n+1,HO}}{\Delta x_i}+{\sigma }_{i}c{\bar{E}}_i^{n+1,HO}={\sigma }_{i}ac(T_i^{n+1}{)}^4$$

$$E_i^{n+1,HO}$$

$${\bar{E}}_i^{n+1,HO}$$

$${\bar{F}}_{i\pm 1/2}^{n+1,HO}$$

$$\frac{{\bar{E}}_i^{n+1,HO}-{\bar{E}}_i^{n,HO}}{\Delta t_n}+\frac{{\bar{F}}_{i+1/2}^{n+1,HO}-{\bar{F}}_{i-1/2}^{n+1,HO}}{\Delta x_i}+{\sigma }_{i}c{\bar{E}}_i^{n+1,HO}={\sigma }_{i}ac(T_i^{n+1}{)}^4+{\mathcal{R}}_i^{n+1},$$

$$\frac{{\bar{F}}_{i+1/2}^{n+1,HO,+}-{\bar{F}}_{i+1/2}^{n,HO,+}}{c\Delta t_n}+\frac{c}{6}\frac{{\bar{E}}_{i+1}^{n+1,HO}-{\bar{E}}_i^{n+1,HO}}{\Delta x_{i+1/2}}+{\sigma }_{i+1/2}{\bar{F}}_{i+1/2}^{n+1,HO,+}$$

$$\frac{{\bar{F}}_{i+1/2}^{n+1,HO,-}-{\bar{F}}_{i+1/2}^{n,HO,-}}{c\Delta t_n}-\frac{c}{6}\frac{{\bar{E}}_{i+1}^{n+1,HO}-{\bar{E}}_i^{n+1,HO}}{\Delta x_{i+1/2}}+{\sigma }_{i+1/2}{\bar{F}}_{i+1/2}^{n+1,HO,-}$$

$${\bar{F}}_{i+1/2}^{n+1,HO,+}$$

$${\bar{F}}_{i+1/2}^{n+1,HO,-}$$

$${\bar{F}}_{i+1/2}^{n+1,HO}$$

$$\frac{{\bar{F}}_{i+1/2}^{n+1,HO}-{\bar{F}}_{i+1/2}^{n,HO}}{c\Delta t_n}$$

$$={\gamma }_{i+1/2}^{n+1,+}c{\bar{E}}_i^{HO,n+1}-{\gamma }_{i+1/2}^{n+1,-}c{\bar{E}}_{i+1}^{HO,n+1}.$$

$$\frac{{\bar{E}}_i^{n+1,LO}-{\bar{E}}_i^{n,LO}}{\Delta t_n}$$

$$\frac{{\bar{F}}_{i+1/2}^{n+1,LO}-{\bar{F}}_{i+1/2}^{n,LO}}{c\Delta t_n}$$

$$={\gamma }_{i+1/2}^{n+1,+}c{\bar{E}}_i^{n+1,LO}-{\gamma }_{i+1/2}^{n+1,+}c{\bar{E}}_{i+1}^{n+1,LO},$$

$$\frac{{\bar{F}}_{i-1/2}^{n+1,LO}-{\bar{F}}_{i-1/2}^{n,LO}}{c\Delta t_n}$$

$$=-{\gamma }_{i-1/2}^{n+1,-}c{\bar{E}}_i^{n+1,LO}+{\gamma }_{i-1/2}^{n+1,-}c{\bar{E}}_{i-1}^{n+1,LO},$$

$$\rho c_v\frac{T_i^{n+1}-T^n}{\Delta t_n}$$

$$\frac{{\bar{E}}_i^{n+1}-{\bar{E}}_i^n}{\Delta t_n}+\rho c_v\frac{T_i^{n+1}-T^n}{\Delta t_n}+\frac{{\bar{F}}_{i+1/2}^{n+1}-{\bar{F}}_{i-1/2}^{n+1}}{\Delta x_i}-{\mathcal{R}}_i^{n+1}$$

$$=$$

$$E_i^{n+1}$$

$$=\frac{1}{c\Delta x}\sum _{p\in i}{w}_p(t_{n+1}),$$

$${\bar{E}}_i^{n+1}$$

$$=\frac{1}{c\Delta x_i\Delta t_n}{\int }_{x_{i-1/2}}^{x_{i+1/2}}{\int }_{-1}^1{\int }_{t_n}^{t_{n+1}}\sum _p^{N_p}{w}_p(t)\delta (x-x_p(t))\delta (\mu -{\mu }_p(t))dtd\mu dx,$$

$$=\frac{1}{c\Delta x_i\Delta t_n}{\int }_{x_{i-1/2}}^{x_{i+1/2}}{\int }_{t_n}^{t_{n+1}}\sum _p^{N_p}{w}_p(t)\delta (x-x_p(t))dtdx.$$

$${\int }_0^{\Delta t}{\sigma }_{i}c{w}_p(t)dt$$

$$w_p(t)$$

Full text

1 Introduction/motivation

Modeling of high-energy density physics phenomena has important applications in astrophysics and inertial confinement fusion [5]. In these applications, thermal X-ray radiation becomes the dominant energy transfer mechanism. Radiative transfer is difficult to model due to the stiff nonlinear interactions between radiation and the host material, which can be described by absorption-reemission physics. Thermal radiative transfer (TRT) without physical scattering can be described by the following equations,

[TABLE]

where $I(\mathbf{r},\boldsymbol{\Omega},\nu,t)$ is the specific intensity of radiation at location $\mathbf{r}$, traveling along direction $\boldsymbol{\Omega}$, with speed $c$ and frequency $\nu$ at time $t$, and $E(\mathbf{r},t)=\frac{1}{c}\int_{4\pi}\int_{0}^{\infty}Id\nu d\Omega$ is the (gray) radiation energy density. Here, $\sigma,~{}\rho,$ and $c_{v}$ are the opacity, density and specific heat of host media, and $B(\nu,T)$ is the Planck function for black body radiation at frequency $\nu$ and material temperature $T$,

[TABLE]

The TRT system is difficult to solve numerically due to its large dimensionality and stiff nonlinear coupling between radiation and material through absorption and reemission. Over the last several decades, the development of efficient numerical algorithms has been an active area of research.

A popular mitigation strategy for the stiff nonlinear coupling is linearization of the absorption-reemission physics, such as the Fleck-Cummings linearization in the Implicit Monte Carlo (IMC) method [12]. Another popular discretization is the discrete ordinate ($S_{N}$) method. The advantage of $S_{N}$ over IMC is that it is able to obtain the asymptotic diffusion limit [18, 20] in a relatively straightforward manner [2, 1].

Another class of discretization techniques, which has gained popularity in recent years, is the method of characteristics (MOC) [11, 3, 15]. MOC solves the transport equation analytically along the particle trajectory (characteristics). MOC is widely used for solving the neutron transport equation, and it is especially effective when accelerated by the coarse-mesh finite-difference (CMFD) method [27]. A similar concept is employed in the MULTI2D radiation-hydrodynamics code [26]. However, the MULTI2D utilizes the ”short” characteristic method, which averages the angular intensity at each surface of the mesh. Recently, an asymptotic preserving, short characteristic method was proposed [28, 29]. The so-called asymptotic preserving unified gas kinetics scheme (AP-UGKS) utilizes nonlinear elimination and the integral form of the transport equation, enabling a time-implicit discretization of the collisional source term without operator splitting between advection and collision. However, because their formulation assumes the cell-face averaged radiative flux being the function of the neighboring cells, the method is limited to using an explicit CFL time-step size. A space-time characteristics method for solving the neutral-particle transport equations has been proposed in Ref. [21]. In their work, the temporal variable is treated as another dimension in order to perform “$n+1$” dimensional ray-tracing.

The method we propose in this paper is called ”deterministic particle”(DP) method. The term “deterministic particle” has been introduced in several previous studies [7, 22, 10, 9]. In these DP methods, the transport equation is solved via operator-splitting. The first step treats the advection (streaming) term analytically via the particle equation of the motion, and subsequently the collision physics are solved. Thus, it introduces operator-splitting error.

Our approach is similar to AP-UGKS and the space-time characteristics method in that it avoids operator splitting and analytically integrates temporal evolution of particle weights. However, our approach combines the DP method with a High-Order, Low-Order (HOLO) algorithm [24, 6] and the formulation is particularly suitable for multifrequency TRT problems. A well-informed reemission source is obtained from the discretely-consistent gray LO system. Because of the presence of the consistent, nonlinear LO system, the contribution from the collision physics (e.g., absorption-reemission) can be analytically integrated without imposing an explicit time-step constraint. The combination of the HOLO algorithm and the DP method makes the overall solution algorithm very efficient.

The rest of this paper is organized as follows. We first show the particle representation of the gray TRT equation in Section 2. Then, our overall TRT solution strategy, based on the HOLO method, is presented. In Section 3, we extend our DP algorithm to the multifrequency problem. Section 4 shows the overall solution algorithm for the developed DP-HOLO method, and highlights the similarities and differences between IMC and MOC. In Section 5, we demonstrate the efficacy of the DP algorithm, and we summarize our findings and future work in Section 6.

2 Methods/Algorithm

For simplicity, we begin with the 1D gray absorption-reemission TRT model,

[TABLE]

We will consider an extension to a frequency-dependent system in the next section. In particle-based methods (e.g., Monte Carlo) the distribution function (i.e., specific angular intensity, $I$) is represented as a collection of particles with weight $w_{p^{\prime}}$ as follows,

[TABLE]

where $x_{p^{\prime}}(t),~{}\mu_{p^{\prime}}(t)$ are the spatial position and direction of the particle $p^{\prime}$ at time $t$, respectively. The evolution equation for the particle weight, $w_{p}$, can be obtained by multiplying Eq. (4) by $\delta(x-x_{p}(t))\delta(\mu-\mu_{p}(t))$ and integrating over phase space to find:

[TABLE]

where, $Q(x_{p},t)=\frac{\sigma acT^{4}(x_{p},t)}{2}$. The formal solution for $w_{p}$ along the characteristics is,

[TABLE]

Note that the particle weight is always positive because the integrand in Eq. (8) is always non-negative. Due to the absence of an acceleration term in Eq. (4), the particle trajectory can be simply expressed as,

[TABLE]

The integral of Eq. (8) can be evaluated analytically when the source is reconstructed with polynomials, as will be the case here. In the following subsections, we describe implementation and algorithmic details of the DP-HOLO method.

2.1 HOLO solution strategy

In order to perform the integral in the particle-weight evolution equation, Eq. (8), the reemission source term must be known a priori. For this, we employ a similar approach to the CMFD-accelerated MOC algorithm for nuclear reactor eigenvalue calculations, but, instead of using CMFD, we use a recently developed, High-Order, Low-Order (HOLO) algorithm [24, 23, 6]. The HOLO algorithm is an iterative, moment accelerated scheme, with the LO system defined by taking the first two moments of Eq. (4) together with the material temperature equation,

[TABLE]

Here, $F(x,t)=\int\mu I(x,\mu,t)d\mu$ is the (gray) radiative flux. Note that we have used the standard $P_{1}$ closure in Eq. (12) instead of the more consistent Eddington factor closure [13]. Solving Eqs. (11) and (12) introduces inconsistencies between HO and LO descriptions, not only at the discrete level, but also in the continuum. Consistency between the HO and the LO system can be regained by adding a consistency term [24] that takes into account the mismatch in truncation errors and the transport effects. Note that there will be no consistency term in Eq. (11) because, as we will show, this equation is satisfied automatically in both HO and LO systems.

2.1.1 Discrete representation of the LO system and HOLO consistency

In the following, we describe how the LO system is discretized consistently with the HO system. We begin with Eq. (11), which is written discretely as,

[TABLE]

where $T^{n+1}_{i}$ is the material temperature at cell $i$, $\Delta t_{n}=t_{n+1}-t_{n}$, $\Delta x_{i}=x_{i+1/2}-x_{i-1/2}$, $\sigma_{i}=\sigma(T_{i})$ and,

[TABLE]

These quantities are tallied from particles during the HO step as described below. Discrete consistency between LO and HO system demands that they have the same local energy balance, Eq. (14). Because the (time-discrete) LO system cannot solve for both $E^{n+1}$ and $\bar{E}^{n+1}$ simultaneously, in order for the LO system to be consistent with the temporal integration of the HO system we rearrange Eq. (14) as:

[TABLE]

where $\mathcal{R}_{i}^{n+1}=\frac{\bar{E}^{n+1,HO}_{i}-\bar{E}^{n,HO}_{i}}{\Delta t_{n}}-\frac{E^{n+1,{HO}}_{i}-E^{n,HO}_{i}}{\Delta t_{n}}$ is the residual source term [8].

Next, we propose the following closure equations for partial radiative fluxes at faces $x_{i+1/2}$ [Eq. (12)]:

[TABLE]

with

[TABLE]

which are also evaluated from particles. Here, we have introduced the consistency term, $\gamma^{n+1,\pm}$, such that the moments of the HO solution always satisfy Eqs.  (19) and (20). Subtracting Eq. (20) from Eq. (19) yields,

[TABLE]

Eq. (24) is the discrete LO closure equation. Finally, the discrete LO system for cell $i$ consists of the following set of equations.

[TABLE]

Discrete consistency can be simply realized because the HO moments satisfy linear Eqs. (18) and (24), given $T^{n+1}_{i}$ and $\gamma^{n+1,\pm}_{i+1/2}$. Thus, upon the convergence of $T^{n+1}$ and $\gamma^{n+1,\pm}_{i+1/2}$, and assuming consistency at the previous time-step, there is a unique solution which satisfies Eqs. (18) and (24) simultaneously. Furthermore, it is also the only solution which satisfies Eqns. (25)-(28).

2.1.2 Energy conservation property

After the nonlinear convergence between the HO and LO systems, summing Eqs.  (25) and (28) gives the following local energy conservation,

[TABLE]

Because we use a single value for the face radiative flux, $\bar{F}^{n+1}_{i\pm 1/2}$, our HOLO scheme is conservative. The overall DP-HOLO algorithm is summarized in Sec. 4.

2.2 Tallying

Eq. (14) requires evaluation of three moments, $E^{n+1},~{}\bar{E}^{n+1},$ and $\bar{F}^{n+1}$, during a DP-HO step. Using the specific angular intensity definition of Eq. (6), the end-of-time step radiation energy density for cell $i$, $E^{n+1}_{i}$, can be evaluated from particles as,

[TABLE]

which is merely the sum of the particle weights located in the cell $i$ at the end of the time step. We note that this interpolation strategy of the particle weight to the mesh corresponds to a “top-hat” function of width equal to the cell size, and is formally first order accurate [4]. The time-averaged radiation energy density at cell $i$ is defined as,

[TABLE]

The above term is computed during the DP step as follows. First, within a time step, $t_{n}<t<t_{n+1}$, we partition the tallying step for each spatial cell, where the material properties are assumed to be constant. Without loss of generality, when a particle passes through the cell $i$ during interval $0<t<\Delta t$, the deposited particle weight in cell $i$ due to particle $p$ is given by,

[TABLE]

Eq. (33) is a more convenient tally choice because it directly represents absorbed energy. In addition, Eq. (33) is required for evaluating energy-weighted opacities (see Sec. 3). We consider next the cases where the reemission source $Q(x,t)$ is either constant or linear in $x$, for which the source integral in Eq. (33) can be performed analytically. For a constant source, the particle weight at time $t$ can be expressed as,

[TABLE]

where $Q_{i}$ is the average reemission source in cell $i$. Note that for IMC with “continuous energy deposition” [12], the particle weight evolves with only the first term in Eq. (34). Consequently, the particle weight quickly vanishes in optically-thick cells, and surface radiative flux tallies become noisy [8]. On the other hand, it is easy to see that the particle weight for the DP method approaches $w_{p}\approx\frac{Q}{\sigma}=\frac{acT^{4}}{2}$ in optically-thick regions, which is the correct equilibrium limit and thus results in accurate surface radiative flux tallies.

Next, multiplying Eq. (34) by $\sigma_{i}c$ and integrating over the time interval $0<t<\Delta t$ yields an explicit expression for the deposited particle weight,

[TABLE]

where $\tau=\sigma_{i}c\Delta t$, and $G(\tau)=1-e^{-\tau}$. The first term in Eq. (35) corresponds to the deposition of the original particle weight. The second term corresponds to the weight gain due to the reemission source, and the third term corresponds to deposition of the reemission source along the particle trajectory.

For a (spatially) linear reemission source within a cell (the one employed here), $x_{i-1/2}<x<x_{i+1/2}$, we may express the spatial variation of the reemission source as,

[TABLE]

Then, along the particle trajectory, $x_{p}(t)=x_{0}+\mu_{p}ct$, we can write the reemission source term as a linear polynomial in time,

[TABLE]

Accordingly, the particle weight at time $t$ can be written as,

[TABLE]

and the weight deposition in the cell becomes,

[TABLE]

This is the formula employed in this study, and is formally second-order accurate.

Lastly, a surface-integrated, outgoing partial radiative flux [Eq. (21)] can be estimated as:

[TABLE]

where $A_{i+1/2}/\mu$ is the effective cross-sectional area.

2.3 Boundary conditions

There are two main types of boundary conditions in radiative transfer simulations: inflow boundary (including vacuum) and reflective boundary. In our DP algorithm, when a particle $p$ escapes from the boundary at time $t$, a new particle $p^{\prime}$ enters the computational domain with the following attributes,

[TABLE]

where $T_{BC}$ is the inflow boundary temperature. As seen from Eq. (43), the only difference between the two boundary conditions is the new particle weight. For either boundary conditions, a new particle “reflects” back into the computational domain. In this way, the particle population remains constant unless particle re-sampling is performed.

The boundary conditions for the LO system are provided in terms of the radiative fluxes. For example, the left boundary radiative flux at $x=0$ can be written as,

[TABLE]

Note that $\bar{F}^{n+1}_{1/2}\equiv\bar{F}^{n+1}(x=0)$. The partial radiative fluxes can be tallied while performing the HO step. Then, we employ the following form of the LO flux boundary condition at the left boundary,

[TABLE]

Similarly for the right flux boundary condition we use,

[TABLE]

2.4 Particle initialization

Particle initialization requires specification of position, angle and weight. We define initial phase-space location of each particle $(x_{p},\mu_{p})$ by uniformly dividing the number of initial particles located in cell $i$, with tensorial product. Number of space $N_{xi}$ and angular points $N_{\mu i}$ ($N_{pi}=N_{xi}N_{\mu i}$) within the cell are input parameters. Then, a particle initial phase-space location in the cell $i$ can be expressed as,

[TABLE]

In addition, each particle has an “effective” phase space volume of $dp=\frac{2\Delta x_{i}}{N_{pi}}$. To determine weight, in this work, we assume that the initial radiation distribution is isotropic in angle, and Planckian in frequency. Then, the initial radiation energy density can be described by the temperature, $T_{0}$ as,

[TABLE]

Using Eq. (6) and the definition of the cell-averaged radiation energy density yields,

[TABLE]

Then, the initial particle weight in the cell $i$ becomes,

[TABLE]

where $N_{pi}$ is the number of particles in cell $i$.

3 Multigroup Extension

Extending the DP algorithm for frequency-dependent problems is straightforward with a multifrequency approximation, in which the group-wise distribution function, $I_{g}$, is an integral value over the frequency range $\nu_{g-1}<\nu<\nu_{g}$,

[TABLE]

With this definition, we can rewrite the TRT equation as,

[TABLE]

where $b_{g}$ is the normalized Planck function for group $g$. Following a similar procedure as in the previous section, we can write the particle weight evolution equation per frequency group as,

[TABLE]

We emphasize here that there is a significant advantage of the DP method vs IMC for multifrequency problems in that the particle trajectory equations, Eqs. (9) and (10), are independent of frequency and thus all frequency weights can be integrated simultaneously. Thus, we can provide phase space resolution with the same number of particles regardless of the number of frequency groups. Although each particle carries additional information, e.g., all group-wise particle weights, $w^{g}$, we can reuse the ray-tracing evaluations for all frequency groups and thus significantly reduce computational effort. This will be demonstrated numerically in Sec. 5.

The gray LO system becomes,

[TABLE]

where the gray opacities are computed as:

[TABLE]

Here, we evaluate $\sigma_{g,i}$ at $T^{n}$, but we use the current radiation energy density and Planckian spectrum for the weights.

4 Overall DP-HOLO Algorithm

Algorithm 1 shows our multigroup DP-HOLO iterative scheme. In this work, we employ a predictor-corrector time-stepping [25]. In our current implementation of the DP-HOLO algorithm, the total number of particles remains constant during a simulation. Any particle that leaves the computational domain immediately reflects back into the domain with an appropriate weight. Optimization and particle-count adaptivity are left to the future work. During each time step, a particle carries the following attributes,

cell id ($i$), 2. 2.

initial/current weight ($w_{0}/w$), 3. 3.

initial/current location ($x_{0}/x$), 4. 4.

initial/current direction ($\mu_{0}/\mu$), 5. 5.

remaining time to the end of time step ($dt$), 6. 6.

phase space volume ($dp$),

We keep the initial particle state in order to allow HO solver restart during HOLO iterations within a time-step. Note that “initial” particle state refers to the state at the beginning of each time step, not the beginning of the simulation.

The LO system is solved iteratively via preconditioned Newton-Krylov method. We utilize the NOX solver package of Trilinos [14]. Because the detail of the LO solver is presented in [23], we only provide a brief summary in here. First, $k$-th Newton-step can be expressed as,

[TABLE]

Here, $\mathcal{F}(\mathbf{U}^{k})$ is the nonlinear residual, evaluated by the current solution $\mathbf{U}^{k}$. $\mathcal{J}^{k}=\frac{\partial\mathcal{F}}{\partial\mathbf{U}}|_{\mathbf{U}^{k}}$ is the Jacobian matrix. The analytical Jacobian matrix is formed at each Newton-step from Eqs. (55)-(58), and solved via a multigrid preconditioned GMRES solver to the relative tolerance of $10^{-3}$. Nonlinear elimination [16, 30] of the material temperature and the radiative flux variables is employed in order to solve the LO system efficiently. Nonlinear elimination is relatively straight forward to implement in the TRT system due to absence of spatial coupling in the nonlinear term (absorption-reemission). The resulting Jacobian matrix for the radiation energy density resembles the Fleck-Cummings linearized equation [12]. We use the nonlinear (relative) tolerance of $10^{-8}$ for checking convergence of the Newton iterations. The computational cost for solving the LO system is essentially independent of the number of particles and time-step size, but depends on the number of the spatial cells. The LO solver typically requires less than 5 Newton iterations and less than 5 GMRES iterations per Newton iteration to converge.

4.1 Algorithm comparison with IMC and MOC

In this subsection, we highlight algorithmic similarities and differences of the DP-HOLO from IMC and MOC. First, all three are particle-based methods, and thus time-integration along particle trajectories can be performed analytically, unlike other grid-based methods such as $S_{N}$.

The first difference between IMC and DP-HOLO is the use of random numbers. As stated previously, our DP-HOLO algorithm does not employ random numbers. In DP-HOLO, each particle is deterministically generated at $t=0$, and it follows a straight trajectory until the end of simulation (or reflects back when the particle reaches the boundary). On the other hand, IMC particles are generated at each time step, based on the volumetric and boundary source terms. The phase-space location of the source particles is randomly sampled. In contrast, source contributions in DP-HOLO are accounted for in the particle weight evolution equation.

Second, IMC linearizes the absorption-reemission term, and resulting in the following linearized TRT equation [12]:

[TABLE]

where $f=\frac{1}{1+\frac{4\sigma ac\Delta tT_{n}^{3}}{\rho c_{v}}}$ is the Fleck factor. IMC introduces an “effective” scattering term via linearization. The particle trajectory must be resampled when the particle undergoes scattering events (in addition to sampling of distance to the scattering/absorption events). It is known that this linearization introduces an additional stability constraint for the time-step size [12, 19]. Instead, DP-HOLO solves nonlinear absorption-reemission directly via the LO system, which allows it to follow the dynamical time scale of the problem (See Sec. 5.1.1).

Lastly, the evolution equation for the IMC particle weight between two events (e.g., absorption/scattering, boundary crossing) only takes into account exponential decay of the weight (the first term of RHS of Eq. (8)). i.e.,

[TABLE]

Thus, IMC particles in optically thick (and cold) regions quickly lose their weight, possibly leading to noisy tallies of end-of-time-step radiation energy density, $E^{n+1}_{i}$.

DP-HOLO is closely related to the space-time-characteristics method (MOC) [21], as both use the same space-time characteristics. However DP-HOLO is specifically designed to ensure discrete consistency and solve multigroup TRT problems. A (temporally and spatially) discretely consistent LO system is key to the successful application of the space-time characteristic method to TRT problems.

5 Numerical Examples

In order to demonstrate effectiveness of the DP-HOLO algorithm, we compare its accuracy and efficiency against the state-of-the-art Implicit Monte Carlo (IMC) method. We note that many accelerations and variance reduction algorithms, as well as novel code optimizations, exist for IMC that we have not considered in our comparisons. Nonetheless, given the several order-of-magnitude performance advantage of DP-HOLO vs. IMC, it is unlikely that such optimizations can alter the fundamental conclusion of this study. For all the results shown in this section, IMC employs a continuous energy deposition (implicit capture), a uniform particle sampling in each spatial cell, and frequency stratified sampling. For DP-HOLO, all the results are obtained with the linear reemission source reconstruction, which has second-order convergence rate. Since we use a point particle (delta-function) to represent distribution functions (formally first-order accurate) and the linear reconstruction for the reemission source (formally second-order accurate), we expect the DP-HOLO algorithm will yield convergence between $O(N_{p})$ and $O(N_{p}^{2})$, depending on the dominant error components.

In the subsequent section, we present convergence plots of IMC and DP-HOLO. Errors are measured by comparing against their own reference solutions obtained from very large number of particles (1280 and 8192 particles per cell for IMC and DP-HOLO, respectively) while fixing the spatial and temporal discretizations. Thus, we are only measuring algorithmic convergence with respect to the number of particles used. We expect the computed error:

[TABLE]

to exhibit a proper measure of convergence rate of both IMC and DP-HOLO. Furthermore, convergence plots are presented with actual CPU time rather than the number of particles because the CPU cost per particle differs significantly between IMC and DP-HOLO, while the CPU time is almost directly proportional to the number of particles in both. Finally, we use 8 particles to discretize the angular variable for all DP-HOLO simulations.

5.1 Gray Marshak Wave

The first example we consider is a gray Marshak wave problem. In this example, we use two different sets of material parameters that correspond to optically thin and thick regions. Table 1 shows the problem definitions. In both examples, the system is initially in equilibrium, $T_{m}=T_{r}=T_{0}$, where $T_{m}$ and $T_{r}=(\frac{E}{a})^{1/4}$ are the material and radiation temperatures, respectively. An incoming radiation source is applied at $x=0$. The opacity is a function of the material temperature as follows,

[TABLE]

Figs. 1-4 show the material and radiation temperature profiles computed with various particle resolutions. The value shown in the legends indicates the sampled number of particles per cell (i.e., the number of reemission source particles per cell in IMC, number of initial particles per cell in DP). Initial time step sizes are $10^{-12}$ and $10^{-11}$ s for Problem 1 and Problem 2, respectively. Time step sizes are increased by a factor of $1\%$ and $5\%$ until they reach the maximum time step size of $10^{-10}$s.

The host medium for problem 1 is optically thick, and the radiation and material temperatures stay in equilibrium throughout the simulation. For the optically thick case, material and radiation temperature are essentially identical. The IMC solutions [Figs. 1(a) and 2(a)] exhibit a noticeable statistical noise in all cases but the one with 1280 particles per cell. On the other hand, the DP-HOLO solutions [Figs. 1(b) and 2(b)] show no significant noise regardless of the number of particles used. Fig. 5 is an efficacy plot for this problem. We use the term “efficacy” to indicate required CPU time for a certain accuracy. As can be seen from Fig. 5, errors in material and radiation temperatures are identical, as expected. We have plotted two sets of data for DP-HOLO. The first one is error vs. HO CPU time (blue and yellow data points), and the other is error vs. total CPU time (black data points). Note that in the total CPU time for IMC and the HO solver CPU time for DP-HOLO are almost proportional to the number of particles used. Thus, it can be seen that the IMC simulations exhibit a typical $O(\sqrt{N_{p}})$ convergence rate, while the convergence rate of DP-HOLO appears to be $O(N_{p})$ (until the error is dominated by other sources such as the temporal and spatial discretizations). For DP-HOLO, the computational cost of the LO solver is about 7 sec, which corresponds to $0.7\%$ to $65.6\%$ of the total CPU time depending on the number of particles as shown in Table 2.

Problem 2 has smaller opacity and incoming radiation temperature, and the material and radiation temperatures are decoupled and differ by about $40$ eV at the wave front. The material temperature profiles do not exhibit statistical noise (unlike problem 1) as seen from Fig. 3(a). However, the IMC radiation temperature profile still shows significant statistical noise, seen in Fig. 4(a), while DP radiation temperature solutions are free from the noise for any number of particles (Fig. 4(b)). Fig. 6 shows the efficacy plot for this problem. Similar to the previous case (Fig. 5), the IMC method converges close to $O(\sqrt{N_{p}})$, while the DP-HOLO method converges between $O(N_{p})$ and $O(N_{p}^{2})$. Data in this figure is arranged in the same way as in Fig. 6. The computational cost of the LO solver is about 10.5 sec, which corresponds to 1.3 % to 76.4% of the total CPU time as shown in Table 3. For both optically thin and thick cases, the DP-HOLO method delivers efficacies many orders of magnitude larger than IMC.

5.1.1 Stability with large time-step

It is well known that the Fleck-Cummings linearization used in IMC can produce nonphysical material temperature solutions when the time step becomes sufficiently large [12, 19]. This unphysical phenomenon is often referred to as the “violation of the maximum principle.” When the maximum principle is violated, the material temperature becomes higher than the boundary (incoming) temperature, and stagnates the radiation wave. The only mitigation in the IMC method is to use a smaller time-step size. On the other hand, because of the nonlinearly consistent LO system, our DP-HOLO algorithm does not violate the maximum principle.

In order to demonstrate this effect, we run the gray Marshak wave problem with time steps of various sizes.

Fig. 7 shows the material temperature profile with the different time step sizes. It is clear that the IMC solution with large time step (e.g., $\Delta t>10^{-10}$s) exhibits unphysical temperature spikes and a slow-moving wave front. On the other hand, the DP solution is virtually indistinguishable for all time step sizes.

5.2 Larsen MG

The next example we consider is a multifrequency Marshak wave problem (a.k.a Larsen’s problem [17]). This problem consists of three regions; region 1: $0<x<2.0$ cm, region 2: $2.0<x<3.0$ cm, and region 3: $3.0<x<4.0$ cm. We use spatial mesh size of $\Delta x=0.05$ cm. Each region has the following frequency-dependent opacity,

[TABLE]

where $\alpha$ is set to $10^{9}$, $10^{12}$, and $10^{9}$ for the regions 1, 2, and 3, respectively. The density and heat capacity are set to $1.0$ [g/cc] and $5.109\times 10^{11}$ [erg/g-eV]. The system is initially in thermal equilibrium at $1.0$ eV and we apply a $10^{3}$ eV Planckian surface source at x = 0.0 cm to start the transient. Each simulation was run to $t=6\times 10^{−10}$ s with $\Delta t=2\times 10^{-12}$s. Frequency groups are split uniformly in the log-scale between $10^{-2}$ eV and $10^{5}$ eV.

Figs. 8 and 9 show the material and radiation temperature profiles computed with various numbers of particles. Compared to the material temperature profile, the radiation temperature is quite noisy when IMC is used. No noisy solutions are observed with the DP-HOLO results. Fig. 10 shows the efficacy plot. Similarly to the previous examples, the convergence rate for the radiation temperature is approximately $O(\sqrt{N_{p}})$ and between $O(N_{p})$ and $O(N_{p}^{2})$ for IMC and DP-HOLO, respectively. Finally, the computational cost of the LO solver is about 6 sec, which corresponds to 0.7 $\%$ to 55 $\%$ of the total CPU time as shown in Table 5.

5.3 Multigroup Efficacy

Lastly, we present the efficacy of the DP-HOLO algorithm in multifrequency problems. First we simply check the algorithmic scaling with respect to the number of groups with a problem with frequency-independent opacity. Here we use the material property of the optically thin gray Marshak wave problem presented previously (e.g., Eq. (67), $\alpha=10^{6}$).

Fig. 11 shows the material temperature profile of the problem with different numbers of groups, which clearly indicates the solution (and hence the physics) is independent of the group structure. Table 6 shows relative CPU increase as a function of the number of groups. The relative CPU time for the IMC simulations in this problem scales linearly with the number of groups. On the other hand, it can be seen that the total CPU time of the DP-HOLO calculation with 64 groups only increases by a factor of 5.2 compared to the single-frequency group case. Furthermore, the CPU cost of the LO (gray) system is almost independent of the number of groups, and the HO CPU time increases only about a factor of 19 for 64 groups. This increase is significantly smaller than the group number increase of 64.

Next we use the Larsen’s problem configuration. Fig. 12 shows material temperature profiles for different numbers of groups at $t=6\times 10^{-10}s$. It shows that, due to the frequency-dependent opacity, the material temperature profiles change significantly when the frequency variable is not well-resolved. (The solution has converged for 64 groups.) Similar to Table 6, Table 7 shows the relative CPU time increase as a function of groups. In this problem, IMC does not scale linearly with the number of groups. This is mostly due to the fact that a single group simulation of this problem becomes optically very thick, and the IMC particles do not travel far (very small cost per particle), which results in artificially small CPU time for the single-group case. On the other hand, the scaling of the DP-HOLO simulations is almost identical to the previous case. This indicates that the CPU cost of DP-HOLO algorithm is almost independent of the physics regimes.

6 Conclusion

We have proposed a moment-accelerated deterministic particle algorithm for multifrequency thermal radiative transfer problems. The DP-HOLO method is closely related to the method of characteristics employed in nuclear reactor calculations. Unlike Implicit Monte Carlo, the particle weight evolution takes into account the reemission source. An explicit expression of reemission source is obtained from a discretely consistent, moment-based low-order description. With numerical examples, we have demonstrated that the DP-HOLO method results in accurate material and radiation temperatures with a very small number of particles compared to IMC. Furthermore, the DP-HOLO method converges between $O(N_{p})$ and $O(N_{p}^{2})$, as opposed to $O(\sqrt{N_{p}})$ for IMC. As a result, DP-HOLO has efficacies (accuracy per unit of cost) several orders of magnitude larger than IMC. Due to the ability of the DP-HOLO particles to evolve weights for all frequencies simultaneously, the advantage of the DP method over IMC increases with multifrequency problems. Furthermore, the relative CPU time increase of DP-HOLO simulation is independent of the physics regime. Ongoing work includes an extension of the algorithm to two-dimensional problems, as well as curvilinear geometries. Future work will focus on developing a new discretization scheme for the LO system that preserves the asymptotic diffusion limit.

7 Acknowledgment

This work was performed with support of the Laboratory Directed Research and Development program at Los Alamos National Laboratory under US government contract DE-AC52-06NA25396 for Los Alamos National Laboratory, which is operated by Los Alamos National Security, LLC, for the US Department of Energy.

8 Reference