A Monte Carlo Method for 3D Radiative Transfer Equations with Multifractional Singular Kernels Christophe GomezOlivier Pinaud

2025-04-27 0 0 2.97MB 34 页 10玖币
侵权投诉
A Monte Carlo Method for 3D Radiative Transfer Equations with
Multifractional Singular Kernels
Christophe Gomez∗Olivier Pinaud†
June 14, 2023
Abstract
We propose in this work a Monte Carlo method for three dimensional scalar radiative transfer equations
with non-integrable, space-dependent scattering kernels. Such kernels typically account for long-range statis-
tical features, and arise for instance in the context of wave propagation in turbulent atmosphere, geophysics,
and medical imaging in the peaked-forward regime. In contrast to the classical case where the scattering
cross section is integrable, which results in a non-zero mean free time, the latter here vanishes. This creates
numerical difficulties as standard Monte Carlo methods based on a naive regularization exhibit large jump
intensities and an increased computational cost. We propose a method inspired by the finance literature
based on a small jumps - large jumps decomposition, allowing us to treat the small jumps efficiently and
reduce the computational burden. We demonstrate the performance of the approach with numerical simu-
lations and provide a complete error analysis. The multifractional terminology refers to the fact that the
high frequency contribution of the scattering operator is a fractional Laplace-Beltrami operator on the unit
sphere with space-dependent index.
Key words. radiative transfer, singular scattering kernels, Monte Carlo method, wave propagation, random
media, long-range correlations.
1 Introduction
Radiative transfer models have been used for more than a century to describe wave energy propagation through
complex/random media [32, 10], as well as neutron transport [40, 51], heat transfer [54], and are still an active
area of research in astrophysics, geophysics, and optical tomography [39, 43, 44, 45] for instance. In this work,
we propose a new Monte Carlo (MC) method to simulate the following radiative transfer equation (RTE)
(∂tu+ˆ
k· ∇xu=Qu,
u(t= 0, x, ˆ
k) = u0(x, ˆ
k),(t, x, ˆ
k)∈(0,∞)×R3×S2,(1)
where S2denotes the unit sphere in R3, and uis the wave energy density in the context of wave propagation or
a particle distribution function in the context of neutronics. The scattering operator Qhas the standard form
(Qu)(x, ˆ
k) = λ(x)ZS2
Φ(x, |ˆp−ˆ
k|)(u(x, ˆp)−u(x, ˆ
k))σ(dˆp),(2)
for σ(dˆp)the surface measure on S2,Φthe scattering kernel, and λ > 0a function modeling the support of the
scattering process. Regions where λ(x)=0are homogeneous and uundergoes free transport. MC methods
have long be used for the resolution of (1), see e.g. [36, 51]. The originality and difficulty in our work lies in
the fact that we consider situations where the mean free time t0associated with Qvanishes in the scattering
regions, that is 1
t0(x)=λ(x)ZS2
Φ(x, |ˆ
k−ˆp|)σ(dˆp)=+∞,where λ(x)>0,(3)
and as a consequence the standard MC representations of udo not apply. Such a scenario arises for instance in
the context of highly peaked-forward light scattering in biological tissues and in turbulent atmosphere, or more
∗Aix Marseille Univ, CNRS, I2M, Marseille, France, christophe.gomez@univ-amu.fr
†Department of Mathematics, Colorado State University, Fort Collins CO, pinaud@math.colostate.edu
1
arXiv:2210.04956v2 [math.NA] 13 Jun 2023
generally in the context of wave propagation in random media with long-range correlations that we describe
below. In this paper we write Φas
Φ(x, |ˆp−ˆ
k|) := a(|ˆp−ˆ
k|)
|ˆp−ˆ
k|2+α(x)=1
21+α(x)/2ρ(x, ˆ
k·ˆp),with ρ(x, s) := ap2(1 −s)
(1 −s)1+α(x)/2s∈[−1,1).(4)
Above, α:R3−→ [0,2) accounts for the slow variations of scattering across the ambient space, and ais a smooth
bounded function characterizing some statistical properties of the medium and such that a(0) >0. Practical
examples are given further. A direct calculation shows that (3) holds when α∈[0,2). Also, the integral in (2)
has to be understood in the principal value sense when α∈[1,2), see [23]. The multifractional terminology
that we use is motivated by the fact that the unbounded operator Qcan be expressed as a (multi)-fractional
Laplace-Beltrami operator (−∆S2)α(x)/2on the unit sphere up to a bounded operator w.r.t. the ˆ
kvariable
[22, 23].
We would like to emphasize that we focus in this work on kernels of the form (4) for simplicity of the
exposition, and that our method applies, after proper decomposition (see [23]), to more general kernels that
behave like (4) at the singularity.
The RTE can be derived from high frequency wave propagation in random media, see e.g. [49]. In such a
context, the velocity field c(x)has the form
1
c2(x)=1
c2
01 + √η V0x, x
η x∈R3, η ≪1,
where c0is the background velocity (that we set to one in the sequel for simplicity), V0is a mean zero random
field modeling fluctuations around the background, and ηis the correlation length of the random medium,
assumed to be small after proper rescaling. The first variable in V0represents the slow variations of the random
perturbations, while the second one corresponds to their high frequency oscillations. The latter are responsible
for the strong interaction between the wave and the medium over sufficient distances. The scattering kernel Φ
is related to the correlation function of V0, and assuming V0is stationary (in the statistical sense) with respect
to the fast variable, a kernel of the form (4) can be obtained from random fields such that
E[V0(x, x′)V0(y, y′)] = pλ(x)λ(y)ZR3
a(|p|)
|p|1+ α(x)+α(y)
2
eip·(x′−y′)dp, (5)
with αranging from 0to 2. Denoting by R(x)the expectation in (5) with y=x,y′=x′+x/η, one can
show that Rbehaves like |x|α(x)−2for |x| ≫ 1, and is therefore not integrable. This is how random fields
with long-range correlations are defined, as opposed to random fields with short-range correlations that exhibit
an integrable correlation function. This approach is of practical interest in biomedical imaging as media with
long-range correlations are able to reproduce experimentally observed power-law attenuations associated with
effective fractional wave equations [20, 25]. The value of the exponents is related to the rate of decay of the
correlation function R, and depends on the nature of the imaged tissues as reported in [14, 26, 27]. Variations
of this exponent can then be used for diagnosis purposes [38, 47].
In Figure 1, we provide examples of such 2D random fields. The top-left picture represents a random
medium with short-range correlations (with a standard Gaussian covariance kernel), while the top-right picture
illustrates a random medium with long-range correlations with α≡1. Because of the singularity at p= 0, one
can observe significantly larger statistical patterns than in the short-range case. In the bottom two pictures,
we highlight the roles of λand α:λcharacterizes scattering regions, and αdefines the correlation structure.
In the inner circle of the bottom-left picture we have α≡0.1, which tends to create shorter range fluctuations
than in the outside where α≡1. In the bottom-right picture, we have a three-layer model for αin which the
inner band exhibits smaller statistical patterns than the outer ones. This type of model is used for modeling
non-Kolmogorov atmospheric turbulences, while standard atmospheric turbulence is modeled with the so-called
Kolmogorov power spectrum
Φ(|k|)∝a(|k|)
|k|11/3,
for |k|in the inertial range of turbulence. This corresponds to the case α= 5/3. This case is not always valid
in experiments as reported in [4, 52, 55], and the statistics of atmospheric turbulence have been shown to vary
with altitude. Models have been derived for instance (see [35] for a review) by considering three ranges (0-2km,
2-8km, and above 8km) with distinct power laws (see Figure 16 for an illustration).
2
Figure 1: Realizations of Gaussian random fields. The upper-left picture represents a field with short-range
correlations, while the upper-right picture depicts a field with long-range correlations with λ≡1and α≡1. In
the lower-left picture, we have λ=1{|x1|<15}, and α= 0.1·1{|x|≤10}+ 1 ·1{10<|x|}. In the lower-right, we have
λ≡1and α(x2)=5/3·1{x2≤2}+ 0.5·1{2<x2≤8}+ 1.9·1{8<x2}.
In the context of biological tissues, the following the Gegenbauer scattering kernel ρGand Henyey-Greenstein
(HG) kernel ρHG are commonly used in the peaked-forward regime [29, 48]:
ρG(x, s) := α g (1 + g2−2g s)−1−α/2
2π((1 −g)−α−(1 + g)−α), ρHG(x, s) := 1
4π
1−g2
(1 + g2−2g s)3/2.(6)
The parameter g∈(−1,1) is called the anisotropy factor, and ρHG is obtained by setting α≡1in ρG. The case
g= 0 corresponds to isotropic energy transfer over the unit sphere, g < 0to dominant transfer in the backward
direction, and g > 0to forward energy transfer. The peaked forward regime is obtained in the limit g→1, for
which 1
(1 −g)αρG(x, ˆ
k·ˆp)∼
g→1
α
2π(2 −2ˆ
k·ˆp)1+α/2=α
2π|ˆ
k−ˆp|2+α.(7)
The case α≡1for the HG kernel is widely used in photon scattering in biological tissues [13, 21, 31]. A
typical realization of the corresponding random field in 2D as g→1is depicted in the top-right panel of Figure
1.
There exist a variety of methods for the resolution of (1) that handle the singular nature of the HG kernel,
see e.g. [18, 19, 33, 34, 37]. They are based on finite differences type discretizations, projections over appropriate
bases w.r.t. the ˆ
kvariable, and approximations of the kernel. Here we propose an alternative approach to handle
singular scattering kernel (4) that is based on a MC method. The latter are popular choices for the simulation
of the RTE when the kernel is smooth, see e.g. [36, 41, 42, 46, 51], essentially for their adaptability to a wide
range of configurations and their simplicity of implementation. A downside is their slow convergence rate, and
there is a vast literature on variance reduction techniques for acceleration. In this work, we focus on the design
of an efficient MC method and postpone any variance reduction considerations to future works.
Our approach is based on an adaptation of a method proposed by Asmussen-Cohen-Rosiński [3, 11] (ACR)
for the simulation of Lévy processes with infinite jump intensity. It relies on a small jumps/large jumps de-
composition of the corresponding infinitesimal generator. The main idea is to approximate the generator of the
small-jump part, which possesses the infinite intensity due to the singularity of kernel, by a Laplace-Beltrami
operator (with respect to the angular variables) on the unit sphere S2. This requires us to simulate paths of a
jump-diffusion process over the unit sphere. For this purpose, we use the characterization of Brownian motion
on the unit sphere given in [5] based on a standard stochastic differential equation (SDE) in R3that is suitable
for space-dependent kernels. This situation is hence more involved than the 2D case we investigated in [24]
3
where the small jumps part can be approximated by Brownian motion on the unit circle for which analytical
expressions are available. Note that, as shown in [24], neglecting small jumps altogether in order to use standard
MC methods leads to large errors, and reducing those comes at significantly increased computational cost.
Denoting by ˆµ(u)the estimator produced by our MC method for some observable µ(u)built on the solution
uto (1), we provide an error estimate of the form
P|ˆµ(u)−µ(u)|> E1+E2+E3≪1
as a theoretical support of our method. Above, E1,E2, and E3are small terms characterizing the various
approximation errors from the original model: the Laplace-Beltrami (i.e. small jumps) approximation, the
discretization error of the diffusion process over the unit sphere, and the MC error. Note that the method we
propose here applies directly to the stationary version of (1)
ˆ
k· ∇xu− Qu=u0,(x, ˆ
k)∈R3×S2,
with source term u0, through the relation
u(x, ˆ
k) := Z∞
0
u(t, x, ˆ
k)dt.
The paper is organized as follows. In Section 2, we introduce probabilistic representations for (1) and its
approximation based on the ACR method. In Section 3, we describe our MC method, state the main theoretical
result regarding the overall approximation error, and detail the simulation algorithms. Section 4 is dedicated
to the validation of the method using semi-analytical solutions. Numerical illustrations are given in Section 5,
where we investigate the role of the strength αof the singularity, both when constant or space-dependent in the
case of non-Kolmogorov turbulence, and compare with solutions for the HG kernel. Section 6 is devoted to the
proofs of our main results and we recall in an Appendix the stochastic collocation method.
The numerical simulations are performed using the Julia programming language (v1.6.5) on a NVIDIA
Quadro RTX 6000 GPU driven by a 24 Intel Xeon Sliver 2.20GHz CPUs station. The codes have been imple-
mented using the CUDA.jl library [8, 9].
Acknowledgment. OP acknowledges support from NSF grant DMS-2006416.
2 Probabilistic representations and approximation
2.1 Representation for (1)
The starting point is the following standard probabilistic interpretation to (1):
u(t, x, k) = Ex,ˆ
ku0(D(t)):= Eu0(D(t)) |D(0) = (x, ˆ
k),
where D= (X, K)is a Markov process on R3×S2with infinitesimal generator
Lf(x, ˆ
k) := −ˆ
k· ∇xf(x, ˆ
k) + λ(x)
21+α(x)/2ZS2
ρ(x, ˆp·ˆ
k)f(x, ˆp)−f(x, ˆ
k)σ(dˆp).
A path, or a realization, of the Markov process Dis often referred to as a particle trajectory. The Xcomponent
of Drepresents the position of a particle, and the component Kits direction. The generator Lcomprises
two terms, the transport part describing free propagation of the particle, and the scattering operator (often
referred to as the jump part in the probabilistic literature) describing the evolution of its direction. The jump
component exhibits a non-integrable singularity leading to a infinite jump intensity and a vanishing mean free
time as expressed in (3).
Note that when λand αare constant, it is shown in [23] that the solution uis unique and infinitely
differentiable in all variables for t > 0for any square integrable initial condition. When λand αare infinitely
differentiable with bounded derivatives at all orders, this result remains valid and we will assume throughout
this work that uis smooth. The same applies to the function uεdefined further in Proposition 2.1.
In order to adapt the ACR method, we introduce the following small region over which the singularity of
the kernel ρ(in (4)) is not integrable, resulting in an unbounded infinitesimal generator L:
Sε
<=Sε
<(ˆ
k) := {ˆp∈S2: 1 −ˆp·ˆ
k < ε}ε∈(0,1).(8)
4
We can now decompose the jump part of the generator Linto two components
Lf(x, ˆ
k) = −ˆ
k· ∇xf(x, ˆ
k) + Lε
<f(x, ˆ
k) + Lε
>f(x, ˆ
k)
:= −ˆ
k· ∇xf(x, ˆ
k) + λ(x)
21+α(x)/2 ZSε
<
+ZSε
>!ρ(x, ˆp·ˆ
k)f(x, ˆp)−f(x, ˆ
k)σ(dˆp),
where Sε
>= (Sε
<)cis the complementary set of region (8) over the unit sphere. The part of the scattering
operator involving Sε
>(with no singularity) is the infinitesimal generator of a standard jump Markov process.
Regarding Sε
<(with the singularity), the following result justifies the approximation of this singular part by a
Laplace-Beltrami operator ∆S2over the unit sphere S2. We will use the notation r′
ε=p1−(1 −ε)2/(2 −ε)in
what follows, and set in the rest of the paper 0< ε ≤ε0<1and 0≤αm≤α(x)≤αM<2.
Proposition 2.1 Let ube the solution to (1) and uεbe the solution to
∂tuε+ˆ
k· ∇xuε=σ2
ε(x)∆S2uε+λ(x)
21+α(x)/2ZSε
>
ρ(x, ˆp·ˆ
k)(uε(ˆp)−uε(ˆ
k))σ(dˆp),
uε(0, x, ˆ
k) = u0(x, ˆ
k),
(9)
for (t, x, ˆ
k)∈(0,∞)×R3×S2, where
σ2
ε(x) := 21−α(x)a(0)πλ(x)
2−α(x)r′
ε
2−α(x).(10)
Assuming a′(0) = 0, for any T > 0, we have
sup
t∈[0,T ]∥u(t, ·,·)−uε(t, ·,·)∥L2(R3×S2)≤ε2−(αM/2) √2T E(u)(11)
where E(u)is defined in (34).
The proof of Proposition 2.1 is postponed to Section 6.1. The term E(u)is independent of εand depends on
derivatives of uw.r.t. ˆ
kup to order 4. Note that the error is of order ε1−(αM/2)/(2−αM)when a′(0) ̸= 0 yielding
a less accurate approximation than for a′(0) = 0. The difference comes from a truncated expansion along the
sphere curvature providing an extra order in εassuming a′(0) = 0. This later assumption holds throughout the
remaining of the paper. Based on (11), we then devise a MC method for (9) instead of (1). The advantage in
using (9) is the fact that the angular diffusion term σ2
ε(x)∆S2is the generator of a Markov process that can be
easily simulated. Indeed, for Wa standard 3D Brownian motion on R3and ×the cross product in R3, it is
shown in [5] that the process Bsolving the SDE
dB =B×dW −Bdt, B(0) ∈S2,
has generator 1
2∆S2. A simple adaptation then gives the desired diffusion coefficient. Since the error is of order
ε2−(αM/2), it is always smaller than ε, and can be adjusted to obtain a desired accuracy. Note also that σε(x)
increases as α(x)gets to 2, and diffusion on the sphere eventually becomes the dominant dynamics.
2.2 Representation for (9)
We interpret (9) as the forward Kolmogorov equation of an appropriate Markov process, and as a consequence
focus on forward MC methods, see e.g. [36] for terminology. Backward equations are simulated in a similar
manner, and can be combined with forward methods for variance reduction techniques [7, 40, 51].
The Markov process we consider for this approach is defined by
Dε(t) := X
n≥0
1[Tn,Tn+1)(t)ψZn
n(t−Tn)t≥0,(12)
where (in the remaining of the paper we extensively make use of the notation z= (x, ˆ
k)):
1. The flow ψz
n= (Xx
n, Kˆ
k
n)is the unique strong solution to the SDE
(dXx
n(t) = Kˆ
k
n(t)dt
dKˆ
k
n(t) = √2σε(Xx
n(t)) Kˆ
k
n(t)×dWn(t)−2σ2
ε(Xx
n(t)) Kˆ
k
n(t)dt, (13)
where ψz
n(0) = z,×is the cross product in R3,(Wn)nis a sequence of independent standard Brownian
motions on R3, and σεis defined by (10).
5
摘要：

AMonteCarloMethodfor3DRadiativeTransferEquationswithMultifractionalSingularKernelsChristopheGomez∗OlivierPinaud†June14,2023AbstractWeproposeinthisworkaMonteCarlomethodforthreedimensionalscalarradiativetransferequationswithnon-integrable,space-dependentscatteringkernels.Suchkernelstypicallyaccountfor...

展开>> 收起<<
A Monte Carlo Method for 3D Radiative Transfer Equations with Multifractional Singular Kernels Christophe GomezOlivier Pinaud.pdf

共34页,预览5页

还剩页未读， 继续阅读

声明：本站为文档C2C交易模式，即用户上传的文档直接被用户下载，本站只是中间服务平台，本站所有文档下载所得的收益归上传人(含作者)所有。玖贝云文库仅提供信息存储空间，仅对用户上传内容的表现方式做保护处理，对上载内容本身不做任何修改或编辑。若文档所含内容侵犯了您的版权或隐私，请立即通知玖贝云文库，我们立即给予删除！

相关推荐

更多
立即下载
分类：图书资源 价格：10玖币 属性：34 页 大小：2.97MB 格式：PDF 时间：2025-04-27

开通VIP享超值会员特权

  • 多端同步记录
  • 高速下载文档
  • 免费文档工具
  • 分享文档赚钱
  • 每日登录抽奖
  • 优质衍生服务

作者详情

相关内容

更多

热门标签

人际关系 配电装置 动力学连接体 力的合成 高考理综 全宋诗作者索引 公务员考试
/ 34
评分 收藏
分享
立即下载
客服
关注