Transcript Fluctuations in Direct Simulation Monte Carlo Alejandro L. Garcia San Jose State University & Lawrence Berkeley Nat.

Fluctuations in
Direct Simulation
Monte Carlo
Alejandro L. Garcia
San Jose State University &
Lawrence Berkeley Nat. Lab.
Numerics for Kinetic Equations, Oberwolfach, Nov. 2008
Hydrodynamic Fluctuations
DSMC,
MD,
DPD,
PIC,
Etc.
Given particle positions & velocities, measure
hydrodynamic variables (density, temperature, etc.)
Fluctuations: Annoyance or Feature?
Part 1 of 2
FLUCTUATIONS AS
ANNOYANCES
MEMs Flows with DSMC
Fluctuations often an annoyance in MEMs simulations.
Need long, expensive calculations to resolve mean values.
Pressure field
2 mm ~ 400 l
Micro-beam simulations require days of massively parallel
supercomputer time to resolve the ~1 m/s flow field
M. Gallis, D. Rader, J. Torczynski (Sandia Nat. Lab.)
Statistical Error (Fluid Velocity)
Fractional error in fluid velocity
Eu 
u
ux

u x2 / S
ux

1
1
SN Ma
where S is number of samples, Ma is Mach number.
For desired accuracy of Eu = 1% with N = 100 particles/cell
1
1
S

2 2
NMa Eu Ma2
S  102 samples for Ma =1.0
S  108 samples for Ma = 0.001
(Aerospace flow)
(Microscale flow)
N. Hadjiconstantinou, A. Garcia, M. Bazant, and G. He, J. Comp. Phys. 187 274-297 (2003).
Sophisticated DSMC
In sophisticated DSMC the time step and cell size vary
dynamically so now Dt and V are also random variables.
Cautionary Example #1
Suppose we dynamically make the cell sizes such
that the number of particles in a cell is exactly N0
A recent paper (J. Comp. Phys. 227, 8035, September 2008)
implements DSMC with roughly this type of sorting.
DSMC Collision Rate
The average number of collisions in a cell is
M  M TRY PACCEPT
Use this
N ( N  1) vr MAX Dt
vr


vr MAX
2V
Want this
N  vr Dt
N ( N  1) vr Dt


2V
2V
?
If N, V, and Dt are correlated
then equality does not hold.
2
For Poisson
N ( N  1)  N
2
Collisions in Cautionary Example
Since the number of particles in a cell is exactly N0
the average number of collisions is
N 0 ( N 0  1) vr Dt 1
M 
2
V
Two problems:  N 0 ( N 0  1)  N 02  N
1
1


V
V
2
Error in Collision Rate
Quick calculation estimates that number of
collisions will be lower by a factor of 1  N
<N>
32
16
8
4
2
Prediction
1.00
1.00
0.98
0.94
0.75
2
Simulation
1.00
Note: In standard
DSMC the
1.00
collision rate is
correct even for
0.99
< N > less than
0.95
one particle per
cell.
0.77
A Quick “Fix” ?
Since the number of particles in a cell is exactly N0
we might think that instead we should compute
the number of attempted collisions as
M TRY
N 02 vr MAX Dt

2V
so that
N 02 vr Dt 1
M 
2
V
Results for Quick “Fix”
Quick calculation estimates that number of
collisions will be higher by a factor of 1  N
<N>
32
16
8
4
2
Prediction
1.03
1.06
1.12
1.25
1.55
1
Simulation
1.03
1.06
Results are
even worse
1.13
than before!
1.27
1.57
Cautionary Example #2
How should one estimate hydrodynamic quantities,
such as fluid velocity, from particle velocities?
Instantaneous Fluid Velocity
Center-of-mass velocity in a cell C
N
mv
i
J
u
 iC
M
mN
Average particle velocity
1
v
N
N
v
iC
Note that u  v
i
vi
Estimating Mean Fluid Velocity
Mean of instantaneous fluid velocity
N (t j )


1
1  1
u   u (t j )  
vi (t j ) 


S j 1
S j 1  N (t j ) iC

S
S
where S is number of samples
Alternative estimate is cumulative average
v (t )



 N (t )
S
u
j
*
N (t j )
iC
S
j
i
j
j
Are these equivalent? If not, which is correct?
Landau Model for Students
Simplified model for university students:
Genius
Intellect = 3
Not Genius
Intellect = 1
Three Semesters of Teaching
First semester
Second semester
Third semester
Sixteen
students
in three
semesters
Total
value is
2x3+14x1
= 20.
Average = 3
Average = 1
Average = 2
Average Student?
How do you estimate the intellect of the average student?
Average of values for the three semesters:
( 3 + 1 + 2 )/3 = 2
Or
Cumulative average over all students:
(2 x 3 + 14 x 1 )/16 = 20/16 = 1.25
Significant difference when there is a correlation between
class size and quality of students in the class.
DSMC Simulations
Measured fluid velocity
using both definitions.
Expect no flow in x for
closed, steady systems
Temperature profiles
T = 4
T = 2
T system
Equilibrium
x
u
20 sample cells
N = 100 particles per cell
10 mean free paths
Anomalous Fluid Velocity
Mean instantaneous
u
fluid velocity
measurement gives an
anomalous flow in the
closed system.
Using the cumulative
mean, u* , gives the
expected result of zero
fluid velocity.
T = 4
T = 2
Equilibrium
Position
Correlations of Fluctuations
At equilibrium, fluctuations conjugate hydrodynamic
quantities are uncorrelated. For example, density is
uncorrelated with fluid velocity and temperature,
 ( x, t )u( x' , t )  0
 ( x, t ) T ( x' , t )  0
Out of equilibrium, (e.g., gradient of temperature
or shear velocity) correlations appear.
Density-Velocity Correlation
Correlation of density-velocity fluctuations under T
 ( x)u( x' )
DSMC
Theory is Landau fluctuating
hydrodynamics

When density is above
average, fluid velocity
is negative
Position x’
u
COLD
A. Garcia, Phys. Rev. A
34 1454 (1986).
HOT
Relation between Means of Fluid Velocity
From the definitions,


N 2  JN
u


u  u * 1

 u *
2 
2


N
m
N


From correlation of non-equilibrium fluctuations,
 ( x)u( x)   x(L  x)T
This prediction agrees perfectly
with observed bias.
M. Tysanner and A. Garcia,
J. Comp. Phys. 196 173-83
(2004).
 ( x)u( x' )
x = x’
Relation with Mechanical Variables
Fluid velocity, in terms of mass and momentum, is
N
mv
i
J
1
iC
u


M
mN
N
N
v
iC
i
so
J
u 
M
Mean Instantaneous
Wrong
u*
J
M
Mean Cumulative
Right
Measured error in mean
instantaneous temperature
for small and large N.
(N = 8.2 & 132)
Error goes as 1/N but only
appears out of equilibrium
where  ( x, t ) T ( x' , t )  0
Mean Inst. Temperature Error
Instantaneous Temperature
Error about 1 Kelvin
for N = 8.2
Position
Predicted error from density-temperature correlation in
good agreement with observed bias in < T >.
A. Garcia, Comm. App. Math. Comp. Sci. 1 53-78 (2006)
Non-intensive Temperature
A
B
Mean instantaneous temperature has bias that goes
as 1/N, so < T > is not an intensive quantity.
Temperature of cell A = temperature of cell B yet
not equal to temperature of super-cell (A U B)
Summary for Part 1
• Fluctuations in DSMC are an annoyance
when measuring hydrodynamic means.
• Modifications of DSMC intended to
improve it can have unexpected
consequences due to fluctuations.
• Instantaneous values of hydrodynamic
variables are biased due to correlations of
fluctuations out of equilibrium.
Part 2 of 2
FLUCTUATIONS AS
FEATURES
Importance of Fluctuations
The calculation of molecular distributions is
fundamental to kinetic theory and hydrodynamic
fluctuations, such as the variance of energy or the
correlation of mass-momentum, are directly
related to these distributions.
Phenomena that are
sensitive to the molecular
distribution (e.g., Arrhenius
law for chemical reactions)
are equally sensitive to
fluctuations.
Energy
Distributions
Activation
Energy
Fluctuations in DSMC
• Hydrodynamic fluctuations (density, temperature,
etc.) have nothing to do with Monte Carlo aspect
of DSMC.
• Variance of fluctuations in DSMC are exact at
equilibrium (Multinomial + Max.-Boltz.).
• Time-correlations correct (at hydrodynamic scale)
• Non-equilibrium fluctuations are correct (at
hydrodynamic scale)
Brownian Systems
“Adiabatic” Piston
Pressure
Hot
P
i
s
t
o
n
“Brownian Box”
Pressure
Cold
Hot
Cold
Chambers filled with particles at different temperatures, equal pressures.
Walls are perfectly elastic yet gases come to common temperature.
Heat conducted by Brownian motion of the piston or of the solid box
Adiabatic Piston by DSMC
Initial State: X = L/4, M = 64 m
NL = NR = 320, TR = 3 TL
1.8
0.55
Temperatures
1.6
0.5
0.5
TR(t)
1.4
0.45
(x-piston)/L
energy per particle
Left
Right
1.2
1
0.35
0.8
TL(t)
0.6
0.4
Piston Position
0.4
0
1
2
3
time
4
Single run
5
time
6
7
8
9
20 run ensemble
0.3
10
5
x 10
0.25
0
1
2
3
time
4
5
time
6
7
8
9
10
5
x 10
Brownian Box by DSMC
Initial State: VL = VR , M = 64 m
NL = NR = 320, TR = 3 TL
0
10
|TL(t)  TL()|
1.8
Left
Right
1.6
TR(t)
Single run
energy per particle
energy per particle
1.4
1.2
1
0.8
eat
-2
10
TL(t)
0.6
0.4
-1
10
0
0.5
1
time
1.5
Ensemble
of 200 runs
2
time
2.5
3
3.5
4
5
x 10
-3
10
0
0.5
1
1.5
time
2
time
2.5
3
3.5
4
5
x 10
Feynman’s Ratchet & Pawl
Carnot*
engine
driven by
fluctuations
Brownian
motors and
nanoscale
machines
* almost
COLD
HOT
Violate
NO. Fluctuations
also lift the pawl,
dropping the
weight back down.
pawl
nd
2
Law of Thermo?
WARM
WARM
Triangula Brownian Motor
Feynman’s complicated mechanical geometry not needed.
Simple geometric asymmetry of Brownian object is enough.
Hot gas
Cold gas
P. Meurs, C. Van den Broeck, and A. Garcia, Physical Review E 70 051109 (2004).
Stochastic Navier-Stokes PDEs
Landau introduced fluctuations into the Navier-Stokes equations
by adding white noise fluxes of stress and heat.
Numerical Schemes
We have investigated numerical schemes for these stochastic PDEs
and find that a three-stage Runge-Kutta scheme works well.
Some tricks are needed to compensate for the smoothing (damping)
of fluctuations by the numerical scheme. For example, interpolation
to cell edges is done as:
Comparison with DSMC
DSMC used in testing the numerical
schemes for stochastic Navier-Stokes
Spatial
correlation
of density
Time
correlation
of density
Shock Drift Benchmark
Standing shock wave
Variance of the shock
position as a function
of time
Shock position
varies as a
random walk so
< x2 >  t
J.B. Bell, A. Garcia, and S. Williams, Physical Review E 76 016708 (2007)
Continuum vs. Particle Methods
Each approach has its own relative advantages
PDE Solvers
DSMC, MD, etc.
• Fast (few variables,
simple relations)
• More Physics
(molecular details)
• Accurate (high order
methods)
• Approximate
(incompressible, inviscid)
• Boundary conditions
• Stable (H-theorem)
Algorithm Refinement mantra is “best of both worlds”
Outline of Algorithm Refinement
Start
Advance Grid
Fill Buffers
Move particles
Overlay to Grid
Reflux
A. Garcia, J. Bell, Wm. Y. Crutchfield and B. Alder, J. Comp. Phys. 154 134-55 (1999)
Fluctuations in DSMC/PDE Hybrid
Correlation of density-momentum
fluctuations in hybrids, compared with
pure DSMC calculation (X-marks)
PDE
DSMC
PDE
Deterministic
PDE hybrid
PDE
PDE
DSMC
DSMC
PDE
PDE
Stochastic
PDE hybrid
Position
S. Williams, J.B. Bell, and A. Garcia, SIAM Multiscale Mod. Sim., 6 1256-1280 (2008).
Position
Fluid Instabilities
Stochastic Navier-Stokes for Rayleigh-Tayor mixing instability.
Horizontal
Slices
Heavy
Vertical
Slices
Light
J.B. Bell, S. Williams, and A. Garcia, submitted (2008).
Time
DSMC Variants for Dense Gases
DSMC variants have been developed for dense
gasses of hard spheres (e.g., CBA, Enskog-DSMC)
and for general potentials (e.g., CUBA).
They are successful in reproducing correct equation
of state (EOS) and transport properties.
However, the hydrodynamic fluctuations are
inconsistent with the EOS (e.g., variance of density
is not consistent with the compressibility).
Stochastic Hard-Sphere Dynamics
vij
vn
vj
vi
DS
D
Properties of SHSD
Stochastic Hard Sphere Dynamics (SHSD) is
equivalent to a fluid with a linear core pair potential.
Pair Correlation
Function ( = 1)
A. Donev, A. Garcia, and B. Alder, Phys. Rev. Lett. 101, 075902 (2008).
Fluctuations are
consistent with the
equation of state
Summary for Part 2
• Fluctuations in DSMC are correct at
hydrodynamic scales and DSMC is useful
for simulations of Brownian motion.
• Stochastic PDE schemes for hydrodynamics
are being developed and DSMC is a useful
benchmark for validating these schemes.
• Dense gas variants of DSMC exist but
difficult to get the fluctuations correct.
DSMC 2009 in Santa Fe
Inn & Spa at Loretto
September 13-19th 2009
RGD 2010 in Asilomar, California
Located on
Monterey Peninsula,
south of
San Francisco