AN OPTIMAL-TRANSPORT FINITE-PARTICLE METHOD
FOR MASS DIFFUSION
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
Abstract.
We formulate a class of velocity-free finite-particle methods
for mass transport problems based on a time-discrete incremental varia-
tional principle that combines entropy and the cost of particle transport,
as measured by the Wasserstein metric. The incremental functional is
further spatially discretized into finite particles, i. e., particles charac-
terized by a fixed spatial profile of finite width, each carrying a fixed
amount of mass. The motion of the particles is then governed by a com-
petition between the cost of transport, that aims to keep the particles
fixed, and entropy maximization, that aims to spread the particles so
as to increase the entropy of the system. We show how the optimal
width of the particles can be determined variationally by minimization
of the governing incremental functional. Using this variational principle,
we derive optimal scaling relations between the width of the particles,
their number and the size of the domain. We also address matters of
implementation including the acceleration of the computation of diffu-
sive forces by exploiting the Gaussian decay of the particle profiles and
by instituting fast nearest-neighbor searches. We demonstrate the ro-
bustness and versatility of the finite-particle method by means of a test
problem concerned with the injection of mass into a sphere. There test
results demonstrate the meshless character of the method in any spatial
dimension, its ability to redistribute mass particles and follow their evo-
lution in time, its ability to satisfy flux boundary conditions for general
domains based solely on a distance function, and its robust convergence
characteristics.
1.
Introduction
Particle methods have attained considerable acceptance in solid and fluid
mechanics, e. g, the Smoothed Particle Hydrodynamics method (SPH) [1, 2],
the Material Point Method [3], the Reproducing Kernel Particle method
[4], the Corrective Smoothed Particle Method [5], the Modified Smoothed
Particle Hydrodynamics method [6], the Optimal-Transportation Meshfree
method (OTM) [7, 8], the dislocation monopole method [9, 10], and oth-
ers. Whereas the intial development of particle methods was often
ad hoc
,
recently there has been increasing awareness of the synergistic connection
between particle methods and optimal transport (cf., e. g., [11, 12, 13]
for background), in as much as the theory of optimal transport provides
a suitable mathematical foundation for formulating particle methods and,
1
arXiv:2305.05315v1 [physics.comp-ph] 9 May 2023
2
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
conversely, particle methods supply a natural discretization of problems of
optimal transport.
The essential structure of
flow
problems is exemplified by Benamou and
Brennier’s reformulation of the compressible Euler flow problem [14, 7] in
terms of two fields, the mass density and the velocity field, which are re-
quired to jointly minimize an action integral subject to a continuity equa-
tion constraint. In particular, the discretization of the problem requires the
formulation of discretization rules for the mass density and for the veloc-
ity field. In the OTM method [7, 8], the mass density is dicretized into
point masses, or
particles
, and the velocity field is discretized by means of
conforming interpolation from nodal values. The convergence properties of
OTM have been established in [15, 16, 17] and applications to solid and
fluid flows, including plasticity and fragmentation, have been presented in
[7, 18, 19, 20, 21], among others.
Fedeli
et al.
[22] applied the OTM method to advection-diffusion prob-
lems. The connection between diffusion problems and optimal transport was
elucidated in the seminal work of Jordan, Kinderlehrer and Otto [23, 24, 25].
In the approach of Fedeli
et al.
[22] the OTM two-field structure is retained,
with the mass density dicretized by means of point masses and the velocity
field, which for diffusion problems is entropic in nature, by means of con-
forming interpolation. A particularly powerful property of optimal transport
methods as they bear on advection-diffusion problems is that they are
ge-
ometrically exact
, in the sense that all advection and volume updates are
computed exactly from the incremental transport mapping. This property
is in contrast with upwinding schemes and Eulerian treatments of the vol-
ume constraint, which are inexact, costly and prone to instabilities. The
geometrical structure of optimal transport (cf., e. g., [26, 13] for relevant
background) is, in fact, a main source of its computational power.
In contrast to the OTM two-field representation of diffusion, so-called
blob methods
work solely with the mass density field, which is discretized
into smoothed mass particles, or
blobs
, of a certain finite width [27, 28]. For
diffusion problems, these finite-particle methods have the appeal of setting
forth single-field numerical schemes, at great gain in simplification. How-
ever, they also raise a number of numerical analysis questions, such as the
consistent variational foundation of the schemes, the enforcement of bound-
ary conditions, the optimal selection of particle thickness and the overall
accuracy and convergence of the schemes.
In the present work, we resort to a time-discretized variational formulation
of mass diffusion, as in [22], but discretize the mass density by means of finite
particles, or blobs, as in [27, 28], in order to formulate a class of velocity-
free finite-particle methods. The immediate benefit of using finite-width
particles, as opposed to Dirac deltas, is that their entropy can be computed
directly from the relative entropy, or Kullback-Leibler, functional (cf., e. g,
[29]). This is in contrast to mass discretization based on point masses,
such as used in [22], for which the entropy functional needs to be carefully
THE FINITE-PARTICLE METHOD
3
rephrased in a weak form. The motion of the finite particles then follows
as the result of a competition between entropy, as given by the Kullback-
Leibler functional, and mobility, as represented by the Wasserstein distance
(cf., e. g., [12]) between consecutive configurations. Thus, entropy works to
disperse the particles uniformly over the domain, whereas mobility works to
hinder their motion.
In Section 3, we show how this competition results in discrete Euler-
Lagrange equations that, in combination with suitable quadrature rules, are
explicit in the positions of the particles. A numerical issue that is central to
finite-particle schemes such as considered here concerns the efficient choice
of particle width. In the present work, we optimize the particle widths
variationally
, i. e., we determine them by minimizing the same incremental
functional that determines the motion of the particles. We show that this
optimality criterion results in particle widths on an intermediate scale be-
tween the mean-free-distance between particles and the size of the domain,
as expected. Another critical issue concerns the formulation of boundary
conditions. We show that flux boundary conditions can be conveniently
formulated in terms of injection/ejection of particles into/from the domain.
We also show that zero-flux conditions can be formulated by means of a
barrier potential based solely on a distance function to the domain, which
greatly facilitates consideration of general domains defined, e. g., by means
of CAD tools.
In Section 4, a number of implementational issues are discussed, includ-
ing the efficient calculation of particle widths and the use of fast searches
to accelerate the calculation of pairwise interaction forces. These issues be-
come critical when the number of particles becomes large and need to be
addressed carefully in order for the calculations to be feasible. In Section 5,
we illustrate the scope and properties of the method by means of a test prob-
lem concerned with the injection of mass into a sphere. The ability to inject
finite particles into the domain and follow the diffusion-driven trajectories
in time, including convergence to a uniform steady-state density after the
injection ends, is quite remarkable. The weak convergence of the scheme is
monitored through the moment of inertia of the particles at steady state.
The test is indicative of robust convergence with respect to the number of
particles.
2.
Problem definition
We consider the advection-diffusion initial-boundary-value problem
∂
t
ρ
+
∇·
(
ρu
) =
κ
∆
ρ
+
ρs,
in Ω
×
[0
,T
]
,
(1a)
ρ
=
g,
on Γ
D
×
[0
,T
]
,
(1b)
(
κ
∇
ρ
−
ρu
)
·
n
=
f,
on Γ
N
×
[0
,T
]
,
(1c)
ρ
(
x,
0) =
ρ
0
(
x
)
,
in Ω
,
(1d)
4
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
where
ρ
is the unknown density or concentration,
ρ
0
is its initial value,
u
is a given advection velocity field,
s
is a source density,
g
is a prescribed
density,
f
is a prescribed outward mass flux,
κ
is a diffusion coefficient, Ω is
a bounded domain in
R
d
, Γ =
∂
Ω is the boundary of Ω,
n
is the outward unit
normal at the boundary, Γ
D
is the Dirichlet boundary and Γ
N
the Neumann
boundary, with Γ
D
∩
Γ
N
=
∅
and Γ
D
∪
Γ
N
= Γ.
Problem (1a) can be equivalently reformulated as the transport problem
∂
t
ρ
+
∇·
(
ρv
) =
ρs,
(2a)
ρv
=
ρu
−
κ
∇
ρ,
(2b)
where
v
is a velocity field that results from the combined effect of advection
and diffusion, eq. (2a) is the continuity equation of Eulerian continuum
mechanics and (2b) is a mobility law combining the effects of advection and
diffusion. The corresponding Lagrangian formulation is
ρ
(
φ
(
x,t
)
,t
) =
ρ
0
(
x
)
det
∇
φ
(
x,t
)
exp
(
∫
t
0
s
(
φ
(
x,τ
)
,τ
)
dτ
)
(3a)
∂
t
φ
(
x,t
) =
v
(
φ
(
x,t
)
,t
)
,
(3b)
where
φ
: Ω
×
[0
,T
]
→
Ω is the transport map. The mass density
(4)
ρ
0
(
x
) exp
(
∫
t
0
s
(
φ
(
x,τ
)
,τ
)
dτ
)
in (3a) gives the initial mass density at material point
x
adjusted for the
mass inserted or removed by the sources. Eq. (3a) states that
ρ
(
x,t
) is the
push-forward of (4) by the transport map
φ
(
·
,t
) whereas eq. (3b) relates the
transport map to the particle velocity.
Eqs. (3a) and (3b) may be thought as jointly defining an evolution for
both the measure
ρdx
and the transport map
φ
. In particular, we note that
this extension to measures requires consideration of the transport map
φ
as
an additional unknown of the problem.
3.
Approximation
We aim to approximate the solutions of (1) variationally through a com-
bination of time discretization followed by spatial discretization of the mass
density into particles. To this end, we resort to an incremental variational
principle that sets forth an optimality criterion for the approximations.
3.1.
Fractional steps.
We begin by discretizing time into a sequence of dis-
crete times
t
0
,t
1
,...,t
k
,t
k
+1
...
and approximating the corresponding mass
densities
ρ
0
,ρ
1
,...,ρ
k
,ρ
k
+1
...
. To this end, we exploit the additive struc-
ture of (1a) to decompose incremental updates by the method of fractional
steps [22]. The advective fractional step is governed by the pure advection
equation
(5)
∂
t
ρ
+
∇·
(
ρu
) = 0
,
THE FINITE-PARTICLE METHOD
5
which follows formally from (2a) by setting
κ
= 0 and
s
= 0 and it be
conveniently solved exactly following the approach described in [22]. This
exact geometrical treatment of advection sets particle methods apart from
upwinding methods and confer them great advantage and strength. The
source fractional step
(6)
∂
t
ρ
=
ρs,
is likewise amenable to an exact treatment, cf. (4). Therefore, in the remain-
der of the section we focus on the pure diffusional fractional step, obtained
by formally dropping the advection and source terms in (2a). Further detail
on the convection step and examples of application may be found in [22].
3.2.
Variational formulation of the diffusional fractional step.
Fol-
lowing [23, 24, 25, 22], we characterize the incremental evolution of the
density by means of the Jordan-Kinderlehrer-Otto (JKO) functional
(7)
F
(
ρ
k
+1
) =
1
2
d
2
W
(
ρ
k
,ρ
k
+1
)
t
k
+1
−
t
k
+
∫
Ω
κρ
k
+1
log
(
ρ
k
+1
ρ
∞
)
dx
→
inf!
,
where
d
W
is the Wasserstein distance (cf., e. g., Appendix A of [22]) and
ρ
∞
is a reference density. For definiteness, we specifically choose
ρ
∞
=
M/
|
Ω
|
,
where
M
is the total mass of the system at steady state and
|
Ω
|
is the volume
of the domain.
It is verified that this formulation indeed supplies a time discretization of
the mobility law (2b) in the limit of zero advection. Thus, a straightforward
calculation (cf., e. g., Appendix A of [22]) gives the Euler-Lagrange equation
(8)
ρ
k
+1
(
x
)
x
−
φ
k
+1
→
k
(
x
)
t
k
+1
−
t
k
=
−
κ
∇
ρ
k
+1
(
x
)
,
where
φ
k
+1
→
k
is the incremental transport map,
φ
k
→
k
+1
=
φ
−
1
k
+1
→
k
and
ρ
k
+1
follows from
φ
k
+1
→
k
through the geometrically exact pushforward operation
(9)
ρ
k
+1
◦
φ
k
→
k
+1
=
ρ
k
/
det
(
∇
φ
k
→
k
+1
)
.
We verify that (8) is indeed a time discretization of (2b). A rigorous analysis
of convergence with respect to the time step may be found in [23, 24].
3.3.
Spatial discretization.
From the standpoint of weak convergence of
measures, it would be natural to approximate
ρ
(
x,t
) by concentrating mass
on a finite collection of points
{
x
p
(
t
)
}
n
p
=1
, i. e., by setting
(10)
ρ
(
x,t
)
∼
n
∑
p
=1
m
p
δ
(
x
−
x
p
(
t
))
,
where
δ
is the Dirac delta and
m
p
is the mass carried by particle
p
. However,
the functional (7) is not defined for trial densities of this type. In [22], this
difficulty is circumvented by replacing the incremental functional (7) with
an alternative weak statement, modeled after (2), in which the mass density
appears linearly and can therefore be approximated by the
ansatz
(10).
6
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
In the present work, following Carrillo
et al.
[27, 28] we instead spread,
or ’fatten’, the particles over a finite width, to be determined. The corre-
sponding approximation is
(11)
ρ
(
x,t
)
∼
n
∑
p
=1
m
p
φ
p
(
x
−
x
p
(
t
)
,t
)
,
where the function
φ
p
(
x
−
x
p
(
t
);
t
) represents the time-dependent profile of
particle
x
p
. The form of the profile
φ
p
(
x
−
x
p
(
t
)
,t
) may be parametrized
within a convenient class of functions and the explicit dependence on
t
then
represents the time-dependence of the parameters. The discrete system
is then obtained by inserting the
ansatz
(11) into the functional (7) and
rendering the resulting discrete functional stationary with respect to all
degrees of freedom.
3.4.
Gaussian particle profile.
We specifically choose a Gaussian particle
profile of the form
(12)
φ
p
(
x
−
x
p
(
t
)
,t
) =
(
β
p
(
t
)
π
)
d/
2
e
−
β
p
(
t
)
‖
x
−
x
p
(
t
)
‖
2
:=
N
(
x
−
x
p
(
t
)
,β
p
(
t
))
,
normalized so that
(13)
∫
R
d
N
(
x
−
x
p
(
t
)
,β
p
(
t
))
dx
= 1
.
We note that
N
(
x
−
x
p
(
t
)
,β
p
(
t
)) has dimensions of an inverse volume, carries
unit mass and has a width set by the parameter
β
p
(
t
).
3.5.
Wasserstein distance.
The Gaussian profile of the particles and the
absence of mass redistribution renders the calculation of the Wasserstein
distance
d
2
W
(
ρ
k
,ρ
k
+1
) in (7) straightforward. Thus, for particles of constant
mass
d
2
W
(
ρ
k
,ρ
k
+1
) simply records the total cost of transportation of all the
particles from their initial to their final positions, namely,
(14)
d
2
W
(
ρ
k
,ρ
k
+1
)
∼
n
∑
p
=1
d
2
W
(
m
p
N
p,k
(
x
)
,m
p
N
p,k
+1
)
,
where we write
N
p,k
(
x
) :=
N
(
x
−
x
p,k
,β
p,k
),
N
p,k
+1
(
x
) :=
N
(
x
−
x
p,k
+1
,β
p,k
+1
),
and neglect the effect of the clipping of the tails of the Gaussians by a finite
domain. Conveniently, the Wasserstein distance between Gaussian measures
is known explicitly [30, 31, 32, 33], whence we obtain
(15)
d
2
W
(
m
p
N
p,k
(
x
)
,m
p
N
p,k
+1
)
∼
(
‖
x
p,k
+1
−
x
p,k
‖
2
+
|
β
−
1
/
2
p,k
+1
−
β
−
1
/
2
p,k
|
2
)
m
p
.
We note that
d
2
W
(
ρ
k
,ρ
k
+1
) quantifies not only the extent of particle displace-
ment during a time step but also the change in the width of the particles, if
any.
THE FINITE-PARTICLE METHOD
7
3.6.
Hermitian quadrature.
The choice of Gaussian profiles also facili-
tates the computation of integrals by numerical quadrature. Thus, for den-
sities of the form (11) and (12) we have
∫
Ω
κρ
k
+1
log
ρ
k
+1
dx
=
n
∑
p
=1
∫
Ω
κ
log
(
ρ
k
+1
(
x
)
ρ
∞
)
×
(
β
p,k
+1
π
)
d/
2
exp
(
−
β
p,k
+1
‖
x
−
x
p,k
+1
‖
2
)
m
p
dx.
(16)
The Gaussian weight in the integral suggests approximation by Hermitian
quadrature [34]. The one-point quadrature rule is particularly simple and
takes the form
(17)
∫
Ω
κρ
k
+1
(
x
) log
(
ρ
k
+1
(
x
)
ρ
∞
)
dx
∼
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
.
Higher accuracy can be achieved by recourse to quadrature rules of higher
order, albeit at increased computational cost and complexity.
3.7.
The discrete problem.
Combining the preceding approximations,
the reduced JKO functional takes the form
F
(
ρ
k
+1
)
∼
n
∑
p
=1
1
2
‖
x
p,k
+1
−
x
p,k
‖
2
+ (
β
−
1
/
2
p,k
+1
−
β
−
1
/
2
p,k
)
2
t
k
+1
−
t
k
m
p
+
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
→
min!
(18)
with
(19)
ρ
k
+1
(
x
) :=
n
∑
p
=1
(
β
p,k
+1
π
)
d/
2
e
−
β
p,k
+1
‖
x
−
x
p,k
+1
‖
2
m
p
.
We note that the incremental variational principle (18) determines both the
updated position of the particles and their optimal widths. The correspond-
ing Euler-Lagrange equations are
(20)
∂F
∂x
r,k
+1
= 0
,
∂F
∂β
r,k
+1
= 0
,
which evaluate to
m
r
x
r,k
+1
−
x
r,k
t
k
+1
−
t
k
+
∂
∂x
r,k
+1
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
= 0
,
(21a)
−
m
r
2
β
−
1
/
2
r,k
+1
−
β
−
1
/
2
r,k
β
3
/
2
r,k
+1
(
t
k
+1
−
t
k
)
+
∂
∂β
r,k
+1
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
= 0
.
(21b)
We note that these equations represent an implicit forward-Euler update for
the positions and widths of the particles.
8
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
3.8.
Boundary conditions.
Neumann, or flux, boundary conditions can
be readily implemented within the particle representation just outlined.
Thus, a prescribed inward flux can be built into the calculations by in-
jecting particles through the boundary of the domain so as to achieve the
desired mass input rate per unit area. A prescribed outward flux can be
implemented either by removing particles crossing the boundary, or by in-
serting particles with negative mass [35], or a combination of both, again so
as to achieve the desired mass input rate per unit area.
A zero flux condition can be implemented by means of a barrier potential,
i. e., a continuous differentiable function
V
(
x
) that is zero on the domain,
including its boundary, and grows steeply away from it. For instance, in
calculations we choose
(22)
V
(
x
) =
C
2
dist
2
(
x,
Ω)
,
where
(23)
dist(
x,
Ω) = min
{‖
x
−
y
‖
:
y
∈
Ω
}
is the distance function to Ω and the constant
C
is to be calibrated. The
incremental problem (18) is then augmented to
(24)
F
(
ρ
k
+1
) +
n
∑
p
=1
V
(
x
p
)
m
p
→
min!
Evidently, the effect of the barrier potential is to apply a restoring diffusional
force
DV
(
x
p
)
m
p
to particles
x
p
that wander out of the domain, causing them
to return to it.
Conveniently, the distance function
d
(
x,
Ω) is a functionality that is pro-
vided by most Computer-Aided Design (CAD) codes, which facilitates the
implementation of the finite-particle method based on CAD models of the
domain of analysis.
A number of extensions of optimal transport and the JKO functional
have been proposed for problems with Dirichlet boundary conditions that
can be applied to particle methods [36, 37, 38, 35]. A possible simple im-
plementation of Dirichlet boundary conditions consists of fixing particles on
the Dirichlet boundary at the desired areal mass density. However, such
extensions are beyond the scope of this paper and will not be addressed
further.
4.
Numerical implementation
In this section, we discuss ways in which the discrete problem (21) can be
simplified, and the solution accelerated, without significantly compromising
accuracy or convergence.
THE FINITE-PARTICLE METHOD
9
4.1.
Computation of diffusive forces.
It is possible to derive a partic-
ular expression for the diffusive forces in (21a) within the accuracy of the
Hermitian quadrature rule. To this end, we write
∂
∂x
r,k
+1
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
=
n
∑
p
=1
κ
[
∂
∂x
r,k
+1
log
(
ρ
k
+1
(
x
)
ρ
∞
)]
x
=
x
p,k
+1
m
p
+
κ
∇
ρ
k
+1
(
x
r,k
+1
)
m
r
,
(25)
with
ρ
k
+1
(
x
) as in (19). Within the accuracy of one-point Hermitian inter-
polation, we have
n
∑
p
=1
κ
[
∂
∂x
r,k
+1
log
(
ρ
k
+1
(
x
)
ρ
∞
)]
x
=
x
p,k
+1
m
p
=
n
∑
p
=1
κ
[
1
ρ
k
+1
(
x
)
∂
∂x
r,k
+1
ρ
k
+1
(
x
)
]
x
=
x
p,k
+1
m
p
∼
∫
Ω
κ
∂
∂x
r,k
+1
ρ
k
+1
(
x
)
dx
=
∂
∂x
r,k
+1
∫
Ω
κρ
k
+1
(
x
)
dx
∼
0
.
(26)
To within this approximation, (25) therefore simplifies to the particularly
simple expression
∂
∂x
r,k
+1
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
∼
κ
∇
ρ
k
+1
(
x
r,k
+1
)
m
r
,
(27)
which we use in all calculations.
4.2.
Uniform particle width.
Considerable simplification is achieved if
the width of the particles is assumed to be uniform, i. e.,
β
r,k
+1
=
β
k
+1
.
Inserting this assumption into (18) and rendering the functional stationary
with respect to
β
k
+1
, (21b) simplifies to the single equation
(28)
−
M
2
β
−
1
/
2
k
+1
−
β
−
1
/
2
k
β
3
/
2
k
+1
(
t
k
+1
−
t
k
)
+
∂
∂β
k
+1
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
= 0
,
where
M
=
∑
n
r
=1
m
r
is the total mass of the system and
(29)
ρ
k
+1
(
x
) :=
n
∑
p
=1
(
β
k
+1
π
)
d/
2
e
−
β
k
+1
‖
x
−
x
p,k
+1
‖
2
m
p
.
10
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
in place of (19). Evaluating the derivative in (28) using (29) yields
∂
∂β
k
+1
n
∑
p
=1
κ
log
(
ρ
k
+1
(
x
p,k
+1
)
ρ
∞
)
m
p
=
Mκ
d
2
1
β
k
+1
−
n
∑
p
=1
κ
(
∑
n
q
=1
‖
x
p,k
+1
−
x
q,k
+1
‖
2
e
−
β
k
+1
‖
x
p,k
+1
−
x
q,k
+1
‖
2
m
q
∑
n
q
=1
e
−
β
k
+1
‖
x
p,k
+1
−
x
q,k
+1
‖
2
m
q
)
m
p
.
(30)
A further simplification is obtained if we assume that the particle widths
are at steady state at all times. Under this assumption, (30) further reduces
to the algebraic equation
M
d
2
1
β
k
+1
=
n
∑
p
=1
(
∑
n
q
=1
‖
x
p,k
+1
−
x
q,k
+1
‖
2
e
−
β
k
+1
‖
x
p,k
+1
−
x
q,k
+1
‖
2
m
q
∑
n
q
=1
e
−
β
k
+1
‖
x
p,k
+1
−
x
q,k
+1
‖
2
m
q
)
m
p
,
(31)
which gives the instantaneous width of the particles as a function of their
positions.
4.3.
Optimal particle width at steady state.
Suppose that, at steady
state, the mass density attains the constant value
ρ
∞
=
M/
|
Ω
|
, where
M
is the total mass and
|
Ω
|
the volume of the domain, assumed to be finite.
Suppose further that this steady state is approximated by a collection of
n
particles of uniform width,
(32)
ρ
(
x
) :=
n
∑
p
=1
(
β
π
)
d/
2
e
−
β
‖
x
−
x
p
‖
2
m
p
,
as in (19). The corresponding relative entropy then follows as
(33)
F
∞
(
β
) =
∫
Ω
κρ
(
x
) log
(
ρ
(
x
)
ρ
∞
)
dx,
to be minimized with respect to
β
. To this end, we rewrite (33) as
(34)
F
∞
(
β
) =
∫
R
d
κρ
(
x
) log
(
ρ
(
x
)
ρ
∞
)
dx
−
∫
R
d
\
Ω
κρ
(
x
) log
(
ρ
(
x
)
ρ
∞
)
dx,
and proceed to estimate each term in turn. We note that in the first term in
(34) the tails of the finite particles are allowed to extend outside the domain,
as tacitly assumed in the quadrature rule (17). The second term in (34) then
corrects for the incurred mass leakoff.
For large
β
the leakoff mass is confined to a narrow layer of thickness
1
/
√
β
in the vicinity of the boundary
∂
Ω of the domain. For sufficiently
regular domains, therefore,
(35)
∫
R
d
\
Ω
κρ
(
x
) log
(
ρ
(
x
)
ρ
∞
)
dx
∼−
κρ
∞
|
∂
Ω
|
β
−
1
/
2
,
THE FINITE-PARTICLE METHOD
11
where
|
∂
Ω
|
is the area of the domain boundary. We note that the effect of
mass leakoff vanishes in the limit of
β
→
+
∞
, as expected. To estimate the
first term of (34), we partition the domain Ω into
n
subdomains of volume
|
Ω
|
/n
, each containing one particle. Again, for sufficiently large
β
, we have
log(
ρ
(
x
))
∼−
β
‖
x
−
x
p
‖
2
, in the area containing particle
x
p
, and we compute
(36)
∫
R
d
κρ
(
x
) log
(
ρ
(
x
)
ρ
∞
)
dx
∼−
κρ
∞
|
Ω
|
β
(
|
Ω
|
n
)
2
/d
.
Inserting these estimates into (34) and minimizing, we obtain
(37)
β
∼
(
|
∂
Ω
|
|
Ω
|
(
d
+2)
/d
)
2
/
3
N
4
/
3
d
,
modulo a multiplicative constant, to be calibrated. In (37), the first factor
is of the form
R
−
2
, where
R
is an effective radius of the domain. The
second factor defines a precise scaling
β
−
1
/
2
∼
N
−
2
/
3
d
for the particle width
with respect to the number of particles. Remarkably, the optimal particle
size takes values on a scale intermediate between the size of the particle
neighborhoods, (
|
Ω
|
/n
)
1
/d
and the size of the domain
|
Ω
|
1
/d
.
4.4.
Local-neighborhood searches.
A na ̈ıve evaluation of the diffusive
forces (27) is an
O
(
n
2
) operation which, for large
n
, becomes the main com-
putational bottleneck. The complexity of the operation can be reduced—and
the calculations accelerated—by noting that, due to the Gaussian decay of
the particle profiles, points
x
p,k
+1
in the sum (19) for the calculation of
∇
ρ
k
+1
(
x
r,k
+1
) at distances much greater than 1
/
√
β
r,k
+1
from the
x
r,k
+1
contribute neglibly to the sum. Therefore, the sum can be restricted to a
local neighborhood of
x
r,k
+1
on the scale of 1
/
√
β
r,k
+1
without appreciable
loss of accuracy.
The rate-limiting step in this approach is the operation of searching for
the local neighborhoods as they evolve during in calculations as a result
of the relative shuffling of the particles. Numerous techniques have been
developed for solving the local-neighborhood problem. One approach makes
use of space partitioning algorithms where the branch and bound method
[39] is applied. Early examples in this class are the
k
-d tree [40, 41] and the
k
-means [42] algorithms, with many variants thereof. A deficiency of exact
methods is that they cannot guarantee logarithmic search time [43, 44].
Conveniently, for many applications approximate nearest-neighbor (ANN)
algorithms [45, 46], which trade off accuracy for efficiency, suffice. More
recent work has focused on proximity graph-based methods [47, 48], which
are well-suited to high-dimensionality problems arising in image recognition,
computational linguistics, product recommendation, and others.
The best choice of search algorithm depends on data-related properties
such as dimension, data-set size, correlations, density distribution, and oth-
ers, and performance metrics such as query time, accuracy, query workload,
building time and memory usage. In this work, we employ a simple
ad
hoc
algorithm based on a regular one-level subdivision of the computational
12
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
Figure 1.
Schematic of the search algorithm with reference to the
problem of an injected sphere. The spherical domain is subdivided in
same-size cubic cells, or boxes, and each particle is assigned to a cell
according to its current position. A cell-connectivity array lists the par-
ticles contained in each cell. For each particle, a second nearest-neighbor
array records the particles contained within a spherical neighborhood of
radius
∼
1
/
√
β
. At every time step, each particle interacts only with the
particles in its neighborhood. Periodically, the nearest-neighbor array
is updated by means of searches restricted to the cells that intersect the
neighborhood of each particle.
domain into cubic cells, or boxes, Fig. 1. For simplicity, we consider the
case of uniform and constant
β
pegged, e. g., to the estimate (37). In this
case, the size of cells is taken as a multiple of 1
/
√
β
, with factors to be
calibrated numerically. Thus, the algorithm requires two connectivity ar-
rays. At initialization, as well as for any new particle entering the domain
as a result of a prescribed flux, the particles are assigned to a specific cell
according to their coordinates. Each particle is outfitted with the connec-
tivity list of interacting particles, and each cell is outfitted with the list of
the particles contained in the cell. In shared memory platforms, the arrays
may be conveniently implemented as C++ STL containers. The algorithm
keeps the cell array up-to-date by transferring the particles between cells at
every time step as needed. By contrast, the particle connectivity lists are
updated only after a prescribed number of time steps. The update is carried
out by searching only the particles contained in the cells intersected by the
spherical neighborhood of the query particle.
THE FINITE-PARTICLE METHOD
13
5.
Example: Mass injection into a sphere
We demonstrate the properties of the finite-particle scheme defined in the
foregoing by means of a test case concerned with the injection of mass into
a spherical domain through a small orifice in its surface.
5.1.
Problem definition.
We consider a spherical domain of radius
R
,
initially empty, which is filled with mass injected through a small orifice of
area
A
at the North pole at constant mass flux
f
during an injection time
interval [0
,T
0
]. Subsequently, the particles diffuse freely, eventually filling
the domain and attaining a uniform steady-state density
ρ
∞
=
M/
|
Ω
|
. The
total injected mass
M
=
fAT
0
is divided into
n
particles of equal mass
m
=
M/n
, each inserted at time intervals of
T
0
/n
to a position below
the North pole offset by a distance
d
from the boundary. The velocity of
injection of the particles is, therefore,
v
=
d/
(
T
0
/n
). The problem data
(dimensionless units) are collected in table 1. To analyze the performance
and the convergence properties of the method, we consider a total diffusion
time
T > T
0
and numbers of particles ranging from a minimum of 125 to
a maximum of 128000. The parameters used in calculations are listed in
Table 1.
R M κ T T
0
v
1 200 1 120 25 1
Table 1.
Parameters used in the sphere-injection test calculations.
The discrete incremental problem (21) is solved by a sequential Runge-
Kutta four-point scheme (RK4) for small systems,
n
≤
8
,
000, and a parallel
implementation of forward-Euler on 1 to 32 cores for larger systems. The
diffusional forces are computed as in (27). For simplicity, the particle width
1
/
√
β
is assumed to be equal for all particles and constant in time, with
value given by the scaling relation (37), calibrated as
β
= 10
n
4
/
9
/R
2
, as
function of the number of particles. The subdivision scheme described in
Section 4.4 is then used with tolerance 0
.
01 to accelerate the calculations.
Owing to the explicit nature of the calculations, the stable time step is
constrained by the diffusional CFL condition ∆
t
∼
∆
x
2
/κ
, where we take
∆
x
∼
(
|
Ω
|
/n
)
1
/d
as an estimate of the distance between particles. Conve-
niently, the stable time step is always smaller than the injection time step.
We additionally enforce the zero-flux boundary condition on the surface of
the sphere by means of the barrier potential (22), with the highest constant
C
that does not decrease the stable time step.
5.2.
Particle flow visualizations.
The evolution of the particle distri-
bution during and after injection is visualized by a sequence of successive
snapshots. Fig. 2 shows the three-dimensional distributions of the centers
of the particles at the end of the injection time
T
0
for four different values
14
A. PANDOLFI
1
, L. STAINIER
2
AND M. ORTIZ
3
of
n
. Every particle is colored according to the number of particles con-
tained in its 1
/
√
β
-neighborhood and its size is related to the actual mass
m
. The distribution of particles is characterized by short-range disorder,
or randomness, and long-range order, or smoothness, and a general trend
towards convergence is evident from the figures. As dictated by optimal
scaling, the particle width increases sublinearly with the number of parti-
cles and the size of the particle neighborhoods also increases concomitantly,
up to a value in the order of 50 particles for
n
= 32
,
000. The value of the
variational principle in elucidating these intricate tradeoffs is remarkable.
(a)
n
=4,000
(b)
n
=8,000
(c)
n
=16,000
(d)
n
=32,000
Figure 2.
Visualization of the 3D distribution of the centers of the
particles at the end of the injection phase
T
0
for different number of
particles. The analysis adopts a constant value of
β
pegged through the
optimal scaling relation to the number of particles
n
. The color of each
particle relates to the number of particles in its neighborhood. The size
of the particles is related to their mass
m
=
M/n
.