Four Decades of Implicit Monte Carlo

ABSTRACT

In 1971, Fleck and Cummings derived a system of equations to enable robust Monte Carlo simulations of time-dependent, thermal radiative transfer problems. Denoted the “Implicit Monte Carlo” (IMC) equations, their solution remains the de facto standard of high-fidelity radiative transfer simulations. Over the course of 44 years, their numerical properties have become better understood, and accuracy enhancements, novel acceleration methods, and variance reduction techniques have been suggested. In this review, we rederive the IMC equations—explicitly highlighting assumptions as they are made—and outfit the equations with a Monte Carlo interpretation. We put the IMC equations in context with other approximate forms of the radiative transfer equations and present a new demonstration of their equivalence to another well-used linearization solved with deterministic transport methods for frequency-independent problems. We discuss physical and numerical limitations of the IMC equations for asymptotically small time steps, stability characteristics and the potential of maximum principle violations for large time steps, and solution behaviors in an asymptotically thick diffusive limit. We provide a new stability analysis for opacities with general monomial dependence on temperature. We consider spatial accuracy limitations of the IMC equations and discussion acceleration and variance reduction techniques.

KEYWORDS:

• Implicit Monte Carlo
• thermal radiative transfer

Appendix

In this Appendix, we discuss modifications to Equation (67a(67a) ) that account for source and sink terms due to a hypothetical weight window treatment to further illustrate how variance reduction mechanisms can be represented in a Monte Carlo transport equation. In this treatment, we suppose that a space-dependent weight window center (or average) is provided as A(x) for x ∈ V, where A(x) ∈ C¹, and 0 < A(x) < ∞. We also are given a so that U(x) = kA(x) is the maximum, upper weight in the weight window and L(x) = A(x)/k is the minimum, lower weight in the weight window. In other words, k is the number of times that a particle will split its weight from w = U to create k particles at weight w = A, and k is the multiplicative weight increase to bring particles from w = L to w = A in a Russian roulette process (assuming that k − 1 particles are killed). The energy-weights are therefore constrained by “boundary conditions” such that L(x) ⩽ w ⩽ U(x).

For simplicity of presentation, we consider Monte Carlo transport processes in which the weight w is constant (no continuous energy deposition or absorption weighting, for example). Therefore, particles stream freely until they undergo collisions or until they intersect with the isosurfaces described by w − L(x) = 0, at which point they undergo Russian roulette, or w − U(x) = 0, at which point they undergo splitting.

We derive the removal term due to particle splitting by considering the net number of Monte Carlo particles removed along a path parameterized by moving from s to s + Δs (along the path ). We presume that we are considering a path along which dU/ds is decreasing so that splitting is possible (if dU/ds > 0, then particles cannot intersect the surface described by w − U(s) = 0) and define ΔU = U(s) − U(s + Δs). Recall from the definition in Equation (61(61) ) that M(s, w) is the number of Monte Carlo particles at s and w per volume, per solid angle, and per energy-weight. Then Recognizing that all of the particles with weights between U(s) and U(s + Δs) will undergo splitting as they stream a distance Δs, we calculate the number of particles removed due to splitting to be (150) Because Δs is assumed small, we may write (151) and, allowing Δs → 0, (152) We therefore identify the net removal rate of particles due to splitting as − (dU/ds)M(s, w = U(s)). Recognizing d/ds as the directional derivative along the particle flight path, the removal term for particle splitting must have magnitude and only occur for weights w = U(x), that is, (153) Another way of deriving this without appealing to number density arguments is to consider an “infinite” splitting probability per unit distance at points , where is defined to be the set of points for which w − U(s) = 0 as ), and to relate the spatial Dirac delta function to the weight-dependent Dirac delta function by (154)

With the removal term so defined, the net rate in the source term due to splitting must match exactly (split particles are only redistributed in energy-weight; their positions and directions are unchanged). In the source term, particles with weights w = U(x) split into k particles with “average” weight A(x) (recall that U = kA), therefore, (155) We have used the definition of to call attention to the detail that may be interpreted as kM particles of weight A or as M particles of weight kA. Incorporating Equations (153(153) ) and (155(155) ) into Equation (67a(67a) ), assuming that Q_m = Q_r = 0 to simplify presentation, we have, (156) Integrating Equation (156(156) ) using ∫^U(x)₀( · )dw, we recover Equation (24a(24a) ), as the splitting source and removal terms cancel: (157)

The source and sink terms describing Russian roulette are entirely comparable; one need only replace L with U and with , as rouletting only occurs along directions of increasing L. We intend to analyze the implications of these mathematical models for various choices of A(x) and k in a future publication.

Notes

1 For instance, in equilibrium, we should find the solution I(ν) = B(ν). Monte Carlo algorithms for TRT treat this problem by sampling from B(ν) directly, not by sampling B(ν)/hνc. Therefore, M(ν) and N(ν) have different spectral distributions.

2 If the term “analog” is defined to mean “a numerical process that is completely analogous to the underlying physics,” then using constant particle weights is not strictly an analog procedure, since Monte Carlo photons used to solve Equation (67(67a) ) have energies that do not depend on their frequency. Regardless, “analog tracking” is often used interchangeably with constant particle weight tracking.

3 Fleck & Cummings (1971) provide details of how to computationally sample a Planckian using only uniform random variates.

4 In Fleck and Cummings (1971), a continuous energy deposition model is used; we discuss this in Section 6.3.

5 This procedure assumes that a constant energy-weight is being used. It is also possible to employ a technique known as absorption weighting or survival biasing, in which a fraction of the particle energy is deposited and the remaining fraction is retained in the radiation.

6 The primary difference between modern implementations of DDMC and IMD is that IMD contains a time-discretization in addition to the spatial discretization (Cleveland and Gentile, 2014).