Transcript lecture5 - Projects at Harvard

Conservation equations
Sauro Succi
Conservation Laws : Hyperbolic equations:
Phi is the (generalized) density and J the corresponding flux;
S is a source (chemistry or other…)
Ubiquituous in all systems with conservation laws
both classical and quantum
Very relevant to all kind of fluids, plasmas, electron transport, you_name_it.
Conservation laws: Gauss theorem
Given a volume V,
enclosed in a piecewise smooth boundary (surface) S,
characterized by the normal n in each point;
Change of “mass” per unit time
= Flux across the enclosing surface
= Flux_in – Flux_out
2
Conservation equations
Hyperbolic Conservation equations:
The flux(current) vector can take many forms, but we shall focus on AD:
HETEROGENEOUS
Advection-Diffusion
Non-linear transport
(self-advecting fluids)
These covers a vast class of HETEROGENEOUS phenomena,
Including Statistical Physics (Fokker-Planck) and Quantum
Mechanics (Schroedinger equation)
Continuity equation
The continuity is a special case with ZERO diffusion and zero source
¶t r + Ñ×(ru) = 0
Looks simpler, but it is not because there is no diffusion,
Sharp features get even sharper, very demanding
on mesh requirements.
Important for all fluids, but especially for near-inviscid ones.
Heterogeneity breaks translational invariance in space (and eventually time).
Non-uniform meshes are usually required to cope with this broken invariance.
In this respect, Finite-Volumes are a more naturalchoice than Finite-Differences
Because they deal with space-averaged quantities instead of point-like ones.
But FD are suitable to insightful and efficient manipulations.
Continuity equation: physics
Number density=Molecules/Volume:
n = N /V
J x (west) = rux DyDz
Mass/time
across surface
mN
¶t r + Ñ×(ru) = 0
r = nm
J x (east) = rux DyDz
Mass change/time = Flux_in – Flux_out
In the limit of zero volume and dt:
Mass density:
The continuity equation
Euler: standing observer
¶t r + Ñ×(ru) = 0
Lagrange: moving observer (“go with the flow”)
Dt (rV ) = 0
Lagrangian (material) derivative
Continuity equation in “ADR” form
Euler: standing observer
Flow compressibility=Source/Sink
of density depending on the flow structure
Mass is ALWAYS conserved!
Poiseuille flow
Incompressible flow in a pipe
ux (x, y) =Uy(1- y); uy (x, y) = 0
Numerics
Looks simple but it is not:
Positivity is a hard constraint!
Numerical Uncertainty Principle
NUP:
At finite resolution it is impossible to
guarantee positivity without a minimum
viscosity
D > Dmin ~ uDx
Why ?:
Wavelengths below dx are missing, but they are
needed to ensure positivity!
Many (eulerian) methods
Euler Space-centered:
unconditionally unstable
Upwind:
overstable, numerical diffusion
Lax-Wendroff:
velocity-dependent numerical diff
Flux-Corrected Transport:
generalized LW for non-uniform velocity profiles
Total Variation Diminishing (TVD) schemes …
Centered Differences
Centered Differences:
n
r n+1
r
j
j
h
u j+1r nj+1 - u j-1r nj-1
= -[
]
2d
Transfer matrix:
Tjk = {+
a j-1
2
,1, -
a j+1
2
}
Troubles with
Stability/Realizability,
One of the two coeff is
necessarily negative!
Upwind
n
r n+1
r
j
j
h
n
r n+1
r
j
j
h
u j r nj - u j-1r nj-1
= -[
]; u j > 0
d
u j+1r nj+1 - u j r nj
= -[
]; u j < 0
d
One sided: first order in BOTH space and time:
Tjk = {a j-1,1- a j , 0}u>0
Tjk = {0,1+ a j ,-a j+1}u<0
Note that a+b+c != 1 for non-uniform profiles…
ok because of compressibility term is a source/sink
Stable, but overdiffusive.
Dnum » ud(1- a )
Lax-Wendroff
Insert a velocity-dependent artificial diffusivity, to lower the
strong numerical diffusivity of upwind.
For homogeneous flows:
U h 2
L = -U¶x + (c +
)¶x
2
2
a=
a
2
(1+ a )
a
b = - (1- a )
2
c = (1- a 2 )
Let us move directly to Flux-Corrected-Transport which generalizes
LW to non-uniform velocity profiles
Flux-Corrected-Transport
Idea: Lagrangian move of a piecewise linear profile over an interval dt,
Interpolate back to the Eulerian grid
J-1
J
J+1
The left side moves by u(j)*dt, while the right side moves by u(j+1)*dt.
Since mass is conserved, the density must decrease/increase according
to the deficit u(j) vs u(j+1).
Between x(j) and x(j+1): INTERPOLATION to get
Rho(J,t+dt) (green oval)
SHASTA
By performing the appropriate algebra:
r
n+1
j
=
n
n
r
r
A
j-1
j
2
-
2
[
d
A
]+ [
2
2
+
r nj+1 - r nj
d
]+ (A- + A+ )r nj
Where:
1/ 2 +a j
A- =
1- (a j-1 - a j )
1/ 2 - a j
A+ =
1+ (a j+1 - a j )
A-2
A-2 + A+2 A+2
Tjk = { , A- + A+ , }
2
2
2
SHarp
And
Smooth
Transport
Algorithm
Homogeneous limit
In the uniform limit U(x)=const:
A- ®1/ 2 + a A+ ®1/ 2 - a
Simple algebra leads to:
r
n+1
j
1 a2 n
= r - [ r - r ]+ ( + )[ r j+1 - 2 r nj + r nj-1 ]
2
8 2
n
j
a
n
j+1
n
j-1
This is centered FD plus a velocity dependent diffusion: Lax-Wendroff
The artificial diffusion is at least 1/8 in mesh units, still a large
number for most applications. Must be improved, how?
Step 2: anti-diffusion
Add an explicit ANTI-diffusion term -1/8*[rho(j+1)-2*rho(j)+rho(j-1)]
But of course this endangers positivity.
To ensure positivity, the antidiffusive fluxes must be LIMITED.
The inspiring criterion is as follows:
The antidiffusive flux canNOT push the value rho(j) above
or below tits neighbors rho(j-1), rho(j+1).
This guarantees that no spurious local extrema (Gibbs!) can
ever be generated by anti-diffusion.
For the detail-thirsty, see
J. Boris and D. Book, J. Comp. Phys. 11,38 (1973),
Several Kcites…
Explicit Anti-diffusion
1 n
A = - [(r j+1 - r nj ) - (r nj - r nj-1 )] º f nj+1/2 - f nj-1/2
8
n
j
This can lead to false oscillations, so
one must secure that No local extrema
arise due to Affusion. Flux Limiters.
n+1
n
r n+1
=
r
+
A
j
j
j
NO!
NO!
j -1
j
j +1
Anti-diffusive fluxes
1 n
A = - [(r j+1 - r nj ) - (r nj - r nj-1 )] º f nj+1/2 - f nj-1/2
8
n
j
1 n
f
= - (r j+1 - r nj ) º AD nj+1/2 ; A = -1/ 8
8
1 n
n
f j-1/2 = - (r j - r nj-1 ) º AD nj-1/2 ; A = -1/ 8
8
n
j+1/2
This is Fick;s law: Flux linearly proportional to the gradient.
It must be turned into a non-linear (ceiling) and non-local law:
Anti-diffusive fluxes: lookaround
D j-1/2 D j+1/2
+
++
¯
++
+
-
+ means positive and ++ means more positive than +
Antidiffusion moves up or down depending on the sign
Of the curvature associated with different combinations:
(+,++), (++,+);(+,-),(-,+),(-,-),(-,--)
The mixed (+,-) and (-,+) should not occurr if one chooses monotonic
Initial Conditions, so there should be only four possibilities.
Step 3: Flux-Limiters: The magic recipe
1
f Cj+1/2 = sgn(D j+1/2 )Max{0, Min{| D j-1/2 |, | D j+1/2 |, D j+3/2 sgn(D j+1/2 )}}
8
1
C
f j-1/2 = sgn(D j-1/2 )Max{0, Min{| D j-3/2 |, | D j-1/2 |, D j+1/2 sgn(D j-1/2 )}}
8
NO!
D j+1/2 = r j+1 - r j
D j-1/2 = r j - r j-1
NO!
j -1
j
j -3/ 2
j -1/ 2
j +1
j +3/ 2
j +1/ 2
Modern developments
Sophisticated forms of FCT which ensure
“Total Variation Diminishing “ (TVD) properties
b
TV[ f ] = ò | f '(x) | dx
a
Obviously very different from
I[ f '] =
b
ò f '(x)dx = f (b) - f (a)
a
TV tells how many slope changes (oscillations) take place in the interval a<x<b.
Whimsical f(x) have large variations even though f(b) may well be equal to f(a)!
By imposing that TV[f] can only decrease in time (like entropy) one
protects against onset of false minima/maxima = Gibbs oscillations
The Fokker-Planck equation
One of the most important equations in Statistical Mechanics,
Chemical-Physics, Quantum Mechanics, but also Social applications
as it deals with Stochastic Processes in general
H. Risken, Springer
The Fokker-Planck equation
One of the most important eqs of math phys: it applies to a broad
variety of stochastic processes in natural (classical and quantum),
life and social sciences
Formally it is a continuity eq with AD fluxes
P(q, t)
Probability density of finding the system around
generalized coordinate q (progress coordinate)
The flux vector consists of a drift (advection)-diffusion term
( Note that: U ~ dq/dt = Progress Rate,
D ~ dq^2/dt)
DNA translocation
0<q<1 is the fraction of translocated DNA
Link to stochastic particle dynamics
(Langevin equation)
Overdamped limit: friction overwhelms inertia:
Advection <-> Force, Acceleration
Diffusion <-> Noise, proportional to the Temperature
U(q) = F(q) / mg
D µ< x > » kT / mg
2
Equilibrium distribution
Steady-state:
This yields:
-1 -F(q)/D
P (q) = Z e
eq
Phi is the generalized potential,
With dimensions L^2/Time=Diffusion
where
Partition Function:
Z(T ) =
+¥
òe
-¥
-f (q)/D
dq
Z contains all the
Thermodynamic
Information!!!
The q-space can be high-dimensional, in which case computing Z
becomes a very non-trivial task (Montecarlo methods)
Maxwell-Boltzmann equilibria
Quadratic potentials
-1 -F(q)/D
P (q) = Z e
eq
Quadratic potential = linear forces:
Gaussian equils:
-1 -kq2 /2D
P =Z e
eq
D / k µ< q2 >
P eq (q)
The particle is gaussian-distributed
around equilibrium position q=0
This is one of the few exceptional cases in which Z can be computed analytically
Phase-Transitions
F(q) = -kq (1- q )
2
2
The attractor q=0 is unstable, the system
jumps erratically from q=-1 to q=+1, double-humped P(q)
Quantum Dynamics in a Bistable Syste
https://www.youtube.com/watch?v=z
Bok
The PDF develops a two-humped structure around q= \pm 1
Complex systems
Corrugated landscapes with competing local
minima (frustration)
Relevant to glasses, protein folding, neural networks…
An elegant and useful map
Fokker-Planck maps 1:1 to Schroedinger equation!
The two equations (1d)
FPE:
¶t P = ¶x [F(x)P + DPx ] = Fx P + FPx + DPxx
SCE: -i¶ty = D q y xx -Vq (x)y xx
Where I call:
F(x) = -U(x)
With a Wick rotation tau=i*t the LHS is the same,
How about the RHS?
It cannot, because SCE has no first order space derivatives.
Can we get rid of it? YES!
FPE versus Schroedinger
Let:
P = gj
Px = gxj + gj x
Where g is an (unknown)
guiding function
Pxx = gxxj + 2gxj x + gj xx
Fx P + FPx + DPxx = Fx gj + F(gxj + gj x )+ D(gxxj + 2gxj x + gj xx )
Collect terms phi, phi_x, phi_xx
(Fx g + Fgx + Dgxx )j + (Fg + 2Dgx )j x + Dgj xx
(Fg + 2Dgx )j x = 0
Namely:
gx / g = -F(x) / 2D
Subscript x
denotes d/dx
FPE versus Schroedinger
The condition:
gives:
where
Remark:
gx / g = -F(x) / 2D
g(x) = Ce
-V (x)/2D
+¥
F(x) º Vx
and
C =1/ ò e-V ( x )/2 D j (x)dx
g (x) º P (x)
2
eq
-¥ (Normalization constant)
Finally, the FPE-SCE mapping is:
Dq = D
Vq (x) = [F 2 (x) / 2D - Fx (x)]+ Dgxx / g
Vq (x) = 1 [(Vx )2 / 2D -Vxx (x)]
2
FPE versus Schroedinger
This is a wonderful result, as it permits to
transfer ADR numerics (and analytics) from
StatMech to Quantum Mechanics, and viceversa!
SUSY Hamiltonians (Fermi/Bose)
V ±SCE (x) = 1 [(Vx )2 / 2D ∓ Vxx (x)]
2
The continuity equation: code
Go to Part1Codes/conti.f and fpe.f
Assignements
1. Solve the 1d continuity equation using
UPWIND and b) SHASTA.
Play with initial conditions (Gaussian, box profile)
and CFL parameters.
2. Derive the flux-limiter expressions
3. Solve the 1d Fokker-Planck equation for
a) Harmonic potential, b) Bistable potential
4. Same in 2D
5. For the brave: same in 3D