Transcript Using Newton Broyden Steps to find SS solutions of the Boltzmann Eq

Acceleration Methods for
Numerical Solution of the
Boltzmann Equation
Husain Al-Mohssen
1
Outline
•
•
•
•
•
•
•
Motivation & Introduction
Problem Statement
Proposed Approach
Important Implementation Details
Examples
Discussion
Future Work
2
Motivation
• Nano-Micro devices have been developed
recently with very small dimensions:
– DLP
– HD read/write head
~ 10  m
~ 0 . 05  m
(Length)
(Gap Length)
• At STP an air molecule travels an average
distance between collisions l  0 . 1 m
• As may be expected the Navier-Stokes (NS)
description of the flow starts to break down as
system length becomes comparable to l
• Accurate engineering models are essential for
the understanding and design of such systems
3
4
Motivation (cnt)
• The Knudsen number is defined as the ratio of
the mean free path to a characteristic
dimension (Kn= l/L). Kn is a measure of the
degree of departure from the NS description
• Kn Regimes:
NS Description Valid
NS Holds inside the domain but slip corrections are
needed at the domain boundaries
Transition Flow
Free molecular Flow
• Recent applications are at low Ma number
5
Introduction
gas
dilute
a
for
A Kinetic Description
 
•
•
A distribution function f ( x , c , t ) is used to describe the gas state, s. t.

 
 
f ( x , c , t ) d c d x is the expected number expected at position x with velocity

c at t.
“Macroscopic” properties are defined as averages over f , for example:


n


3
fdc ;
u 
 c x fdc
3

•
Evolution of f is governed by the Boltzmann Equation
•
Air at STP is satisfies the dilute gas criterion ( n  3  1)
6
Introduction (cnt)
The Boltzmann Equation (BE) in normalized form:
f
  f
 f


c .   a .   Collision
Dt
t
2
x
C
Df
Collision
Integral 

2

Integral

( f ` f 1 `  f f 1 )V  d  d c
2
3
• Follows from the dilute gas assumption
• Valid for all Kn
• 7D(1time+3Space+3Velocity) nonlinear Integrodifferential equation
7
Introduction (cnt)
Numerical Methods of Solving the BE:
• Particle based: DSMC
– Collisionless advection step + collision steps are successively applied.
– Can be shown to simulate BE exactly in the limit of large numbers
[Wagner 1992].
– Chronic sampling problems at low speeds [Hadjiconstantinou et al,
2003].
» Low Ma lmit particularly troublesome
• Approximations of the BE
– Linearized (has many advantages espcially when Ma<<1; still requires
numcerical solution)
eq
– BGK CI Replaced with ( f  f ) / t
• Numerical solutions of the BE
– Recently Baker and Hadjiconstantinou (B&H) proposed a method to
solve the BE at low Ma in a relatively efficient manner.
8
Introduction (cnt)
B&H method of calculating the collision integral:
• Solves the nonlinear BE exactly
• f is written as f  f MB  f D
•
f
MB
f
D
Since f
f
MB
is Maxwell-Boltzmann equilibrium distribution and
is deviation from MB distribution
is not changed by BE, effort is spent on solving
D
•
Even when f D is large the solution is still correct only less
efficient.
• Solution has constant relative noise that is quite small in
contrast to DSMC
B&H solution methods for f:
– Explicit time integration scheme:
• uses time splitting to apply convection step and collision
step separately
• Stability condition limits us to relatively small time steps
– Implicit scheme for finding steady state solutions:
• Scales badly with lower Kn.
– New proposed method for finding SS solutions
9
Problem Statement
We want to find the steady state solution for the first few moments
of the BE (velocity, temp, etc.)
•
•
Consider the x-direction flow velocities in the plot and let us denote u i the
velocity at node i in a certain time

Furthermore, let u (t ) be the vector

T
u ( t )  {u 1 , u 2 ,....., u i ,....... u n }
•
If we define
our system
 

F (u )   u /  t

 
F ( u ss )  0

then we are interested in finding u ss such that for
10
Proposed Solution Methodology
•
•
•

 
We will solve the (in general) nonlinear system of equations F ( u )  0
using Newton’s Method.
In 1-D, Newton’s Method finds successive approximations to F(u)=0 using
F(u) and dF/du=F’(u)
Analogously in multi-dimensions:
F(ui) and F’(ui)
 


1
u i 1  u i  [ J i ] F ( u i )
F(u)
Where the [ J i ] is the Jacobian
matrix of partial derivatives
•
Each iteration of the method will need to evaluate

x
ui+1
ui
 
F (u i )
and the
corresponding [ J i ] to find u i 1 . Since the Jacobian matrix is large and very
•
expensive to compute, a method to approximate new [ J i ] efficiently has
to be found for this approach to be practical
Broyden [Broyden] developed an update method that is very powerful
11
Proposed Solution Methodology (cnt)

The Broyden update formula is a method of updating [ J i ] to [ J i 1 ] such
that:
 
 [ J i 1 ] will be consistent with the new “measured” F ( u i )

[ J i 1 ] will retain as much information as possible from [ J i ] .

Using the Broyden update formula each Newton iteration will only need an
 

evaluation of F ( u i ) to get a new guess of the solution u i 1 and a new [ J i 1 ]

In 1D, Broyden’s method reduces to the Secant Method
12
Simplified Flow Chart of Method
Start
Find
[J ]

Estimate u 0
Integrate BE to find
 
F (u i )
Use Broyden to find [ J i 1 ] from [ J i ]
Find
No
 
and F ( u i )

u i 1
Converged?
Yes
End
13
Important Implementation Details
(for Broyden Portions)


Finding an Initial Jacobian Matrix
 Use continuum solution approximation [ J c ]
 Fairly robust even when [ J c ] is not close to [ J exact ]
Noise
 
o Due to the statistical nature of the method the value of F ( u i ) will
have a noisy component
 


o We can easily show that | x Br  x Ex | N   { F ( x )}

x Ex Is exact solution

x Br Is solution after many Newton-Broyden steps
N  is system characteristic time constant (in steps).
 Less noise is needed for systems with larger time
constants if we want to maintain solution accuracy.
14
1D Graphical Analog
F[u]
 
 { F ( x )}  input noise
 


| x Br  x Ex | N   { F ( x )}


| x Br  x Ex | uncertaint
y in sol
u
15
Important Implementation Details
(BE Portions)

 
Initialization of a proper f ( x , c ) to use with a certain u i :

 
• A finite number of moments at any x is not enough to specify f ( x , c ) at
that position
 
• Need to find a special f ( x , c ) consistent with BE dynamics and

prescribed u i
 
• Use BE to “mature” f ( x , c ) for us by:
1.
2.
3.
4.
•
starting with f MB
 
Integrating BE for a short time to get f ' ( x , c )

Modify f to keep the new shape but give the target u
Repeat 2 and 3 until we converge
Notes:
•
•
Takes 1-4  molecular to converge
 
Has to be done at every evaluation of F ( x i )
1
0.6
0.6
0.8
0.5
0.5
0.4
0.4
0.6
Integrate
BE
0.4
0.2
-1
0.3
1
1
2
3
0.1
-1
0.3
Shift f to
target mean
0.2
1
2
2
3
0.2
0.1
-1
1
3
2
16
3
Flow Chart of Method
Start
Find
[Jc ]

Estimate u 0
Integrate BE
Step1: Equilibrate f
Step2: Sample
Calculation
 
to find F ( u i )
Use Broyden to find
[ J i 1 ]
from
[Ji]
 
and F ( u i )

Find u i 1
No
Converged?
Yes
End
17
Examples
•
•
•
Method successfully calculates 1D flows over the Kn Spectrum (both
pressure driven and shear driven).
Next results are for shear driven flows with a  0.05 (normalized) wall
velocity at different Kn and discretization.
Plot on right shows error bars for different discretization for a kn=0.1
Shear flow. Accuracy of solution is well within expected bounds
U
0.04
Exaggerated
Kn Layer
0.02
100
200
300
400
500
Node #
-0.02
-0.04
Exaggerated Kn
Layer
18
Examples (cnt)
Knudsen Layer
0.0015
kn=0.1
0.001
Broyden Solution
0.0005
Exact layer
20
40
60
80
100
120
-0.0005
-0.001
-0.0015
Convergence History
-2.5
-2.75
-3
512 nodes, kn =0.1
-3.25
-3.5
-3.75
19
2.5
7.5
10
12.5
15
17.5
20
Discussion
•
•
To first order, cost is about time to integrate O(10) iterations. Which
translates to the time to integrate the system about 40  mol (where  mol is
the mean time between collisions)
 
To find solution of accuracy  tg , the allowed error in estimating F ( x i )
(which we will denote  Br ) should be  Br 
•
 tg
N
To increase accuracy, more accurate sampling is required
20
Future Work

Reduce Broyden Integration time:
o Reducing sampling steps by better understanding which
parameters affect the noise level the most.
 
o Refine relaxation method for f * ( x , c )
 able to use method on lower Kn systems

Extend our approach to do 2D/3D grids which would allow more complex
problems with staggered timescales
21
The End
Questions?
22
DSMC Performance Scaling
•
Noise in DSMC is well understood [Hadjiconstantinou et al] and scales as
1
in general
# of Samples
•
It can be shown that :
• Time to find solution by direct integration
 N  Log (
•
 tg
 1
)

 tg




2
Time to find solution by Broden method for similar accuracy
N
 

 tg
•
1




2/3
Log (
1
 tg
)
At large enough N  Broyden method can be significantly faster than direct
integration using the DSMC to solve the BE. This however is only the case
for fairly large N  (of order 104 105)
23
B&H Performance Scaling
•
•
Direct integration cost scales in a similar way to DSMC
Broyden methods performance scales in a more complex manner:
– B&H noise verses cost scaling is more complicated than:
– Noise= fun ( N Samples ,  t ,  x , C . I . parameters ,....) and is
–
–
–
generally nonlinear.
Noise vs. cost scaling becomes similar to DSMC
scaling but only in the limit of very low levels of noise
(& fine meshing). Our numerical experiments indicate
that it is in general much better behaved.
Advantage is even stronger when looking at problems of
engineering accuracy
In our runs Kn~0.1 breakeven point vs. direct integration
24
Plot of Convergence Rates of
Different Methods
• Plot of error for Direct integration, Broyden
and Baker Implicit code. Kn=0.025 # of
nodes 128. (log[Error] vs. log[CI
evaluations])
25
Error of Broyden vs. noise of F
• Show how sig=sig/N_inf in
multidimensions
26
Broyden Step
• Broden formula
• Formula constraints
• Broyden Formula derivation
27
Backup slides+notes
• [[check conv. History 4 high kn and 512]]
• “proper” kndsen layer with 100^3 and
lower noise kn=0.1 and at least 128
nodes. Replace one already in
presentation
• Change Conv. History plto to 512 and
kn0.025 and 30^3 cells
• N_inf vs. Kn for our pb’s to show our rough
break point….
28
DSMC Performance Scaling
(Backup)
Direct Integration Cost:
Broyden Cost:
Slope Sampling Scaling is key:
Sample
2
Ns3
Ns
2
Analysis assumes sampling a small portion of run =>
N
4
3
tg
29
B&H Noise for Different
Paramters(Backup)
For little extra computational Effort you get a dramatic decrease in
measurement error. compare for example pt. A, B and C.
dt.01,g.1Red
dt.1g.1Green dt.1g90Blue
dt.1,.01g1Orange
0.0001
A
Kn=?
0.00001
If only interested in eng.
Accuracy
N_inf=10^-4/sig_sample
B
1.´ 106
1.´ 107
Cost A=Cost B
Cost C=10 Cost A
C
1.´ 108
10
100
1000
30
Distribution Function initilization
(Backup)
• Plot of norm f vs. step [[Possibly for
multiple kn
1
0.8
0.6
[[what kn? What state of F?]]
0.4
0.2
31
1000
2000
3000
4000
Scaling Arguments (Backup)
•
Why is it always O(10)? Well possibly because of this:
•
As per Kelly Newton’s is q-Quadratic and secent is Q-superlinear; Broyden is
somewhere in between.
The other plot is the MMA result using [a] x/nnn + noise
Kelly says eps=K eps^2 not exp[-2t]
•
•
MMA Model Problem in Multi-D
with Noise
-6
-6.2
-6.4
-6.6
-6.8
32
2.5
7.5
10
12.5
15
17.5
20
Can u answer these Questions
• Is it possible that O(10) will increase with
less noise Requrement
• If u reduce Dt sample to decrease noise,
don’t u increase N_inf??!!!
• [[Re-initializing a Run after it reaches its
minimum noise level with less noise as a method
of Confirming convergance or reducing noise
(NB: since we are somehow finding the null
space of the Jacobian aren’t we somehow
garanteed to have a sick matrix when we stall?)]]
33
Can u Explain B&H?
• What is importance sampling? & how is it
applied to CI? Write the appt. version of
CI.
• What is control variate M/C interation?
• How is the finite volume Spliting method
implemented? What are the various
Stability conditions?
34
Integration Stability Codnition
• CI step
• Convection Step
• Implicit step?
35
-7
-6
-5
-4
Log10[Input Noise]
-1
Log10[error]
-2
32
-3
16
512
-4
128
-5
64
36