On Direct and Adjoint Lattice Boltzmann Equations Fran cois Dubois - - PowerPoint PPT Presentation
On Direct and Adjoint Lattice Boltzmann Equations Fran cois Dubois - - PowerPoint PPT Presentation
ICMMES Conference, Hong Kong, July 25-29, 2005 On Direct and Adjoint Lattice Boltzmann Equations Fran cois Dubois Numerical Analysis and Partial Differential Equations Department of Mathematics, Universit e Paris Sud conjoint work with
On Direct and Adjoint Lattice Boltzmann Equations
Scope of the lecture
1) On Taylor and Chapman Enskog expansions for discrete Boltzmann dynamics 2) Inverse methodology for a linear thermal problem 3) Adjoint Lattice Boltzmann Equation
ICMMES Conference, Hong Kong, July 28, 2005
General framework
x : a node of the lattice ∆t : the (small) time step vj : discrete celerity in the lattice, components vα
j , α = 1, 2
then x + vj ∆t is an other node of the lattice f j(x , t) : density of particles having the velocity vj at node x and at discrete time t . m0 = ρ ≡
- j
f j : density of matter mα = qα ≡
- j
vα
j f j : component number α of the momentum
W ≡ ( ρ , q1 , q2 ) : conserved variables
ICMMES Conference, Hong Kong, July 28, 2005
General framework
1 3 6 5 2 4 8 7
Neighbourhood of a vertex x for the D2Q9 LBE model
ICMMES Conference, Hong Kong, July 28, 2005
d’Humi` eres (1992) representation of the discrete dynamics
mk ≡
- j
M k
j f j (k ≥ 3) : other components of the momentum
remark that we have M 0
j = 1 and M α j = vα j
mk
eq and f j eq : equilibrium momenta and velocity distribution
- collision step
mi
∗ ≡ mi ≡ mi eq : conserved momenta (i = 0 , 1, 2) during the collision
mk
∗ = (1 − sk) mk + sk mk eq ,
k ≥ 3 : nonconserved momenta Classical stability condition for the explicit Euler scheme : 0 < sk < 2 f j
∗ ≡
- k
(M −1)j
k mk ∗ :
particle distribution after the collision
- advection step
f j(x , t + ∆t) = f j
∗(x − vj ∆t , t)
ICMMES Conference, Hong Kong, July 28, 2005
Taylor expansion at the order zero
f j(x , t + ∆t) ≡ f j
∗(x − vj ∆t , t)
f j(x , t + ∆t) = f j(x , t) + O(∆t) f j
∗(x − vj ∆t , t) = f j ∗(x , t) + O(∆t)
then mk
∗ =
- j
M k
j f j ∗ = mk + O(∆t)
mk
∗ − mk = O(∆t)
but mk
∗ − mk ≡ −sk (mk − mk eq)
and sk > 0 if k ≥ 3 thus mk = mk
eq + O(∆t)
mk
∗ = mk eq + O(∆t)
coming back in the space of velocity distribution : f j = f j
eq + O(∆t)
f j
∗ = f j eq + O(∆t) .
ICMMES Conference, Hong Kong, July 28, 2005
Taylor expansion at the first order
f j(x , t + ∆t) ≡ f j
∗(x − vj ∆t , t)
f j(x , t + ∆t) = f j(x , t) + ∆t ∂tf j + O(∆t2) f j
∗(x − vj ∆t , t) = f j ∗(x , t) − ∆t vβ j ∂βf j ∗ + O(∆t2)
we take the moment of order k of this identity mk + ∆t ∂tmk + O(∆t2) = mk
∗ − ∆t
- j
M k
j vβ j ∂βf j ∗ + O(∆t2)
we use the previous Taylor expansion at the order zero mk + ∆t ∂tmk
eq + O(∆t2) = mk ∗ − ∆t
- j
M k
j vβ j ∂βf j eq + O(∆t2)
- k = 0 :
∂tρ + ∂β qβ = O(∆t) conservation of mass introduce the tensor F α β ≡
- j
vα
j vβ j f j eq
- k = α :
∂tqα + ∂β F α β = O(∆t) conservation of impulsion
ICMMES Conference, Hong Kong, July 28, 2005
Taylor expansion at the first order
mk + ∆t ∂tmk
eq + O(∆t2) = mk ∗ − ∆t
- j
M k
j vβ j ∂βf j eq + O(∆t2)
but mk − mk
∗ ≡ sk (mk − mk eq)
if k ≥ 3 introduce θk = ∂tmk
eq +
- j
M k
j vβ j ∂βf j eq ≡
- j
M k
j ( ∂tf j eq + vβ j ∂βf j eq )
for i = 0, 1, 2, θi = O(∆t) : Euler equations of gas dynamics then for k ≥ 3 : mk = mk
eq − ∆t
sk θk + O(∆t2) mk
∗ = mk eq −
1 sk − 1
- ∆t θk + O(∆t2)
∂βf j
∗ = ∂βf j eq − ∆t
- k≥3
1 sk − 1
- (M −1)j
k ∂βθk + O(∆t2)
ICMMES Conference, Hong Kong, July 28, 2005
Taylor expansion at the second order
f j(x , t + ∆t) ≡ f j
∗(x − vj ∆t , t)
f j(x , t + ∆t) = f j(x , t) + ∆t ∂tf j + 1
2 ∆t2 ∂2 ttf j + O(∆t3)
f j
∗(x − vj ∆t , t) = f j ∗(x , t) − ∆t vβ j ∂βf j ∗ + 1 2 ∆t2 vβ j vγ j ∂2 βγf j ∗ + O(∆t3)
we take the moment of order i (0 ≤ i ≤ 2) of this identity mi + ∆t ∂tmi + 1 2 ∆t2 ∂2
ttmi + O(∆t3) =
= mi
∗ − ∆t
- j
M i
j vβ j ∂βf j ∗ + 1
2 ∆t2
j
M i
j vβ j vγ j ∂2 βγf j ∗ + O(∆t3)
we use the microscopic conservation mi
∗ = mi
and the previous Taylor expansion at the order one : ∂tmi + 1 2 ∆t ∂2
ttmi + O(∆t2) = −
- j
M i
j vβ j ∂βf j eq +
+ ∆t
- j, k≥3
M i
j vβ j
1 sk −1
- (M −1)j
k ∂βθk + 1
2 ∆t
- j
M i
j vβ j vγ j ∂2 βγf j eq + O(∆t2)
ICMMES Conference, Hong Kong, July 28, 2005
Taylor expansion at the second order
∂tmi +
- j
M i
j vβ j ∂βf j eq = ∆t
- j, k≥3
M i
j vβ j
1 sk − 1
- (M −1)j
k ∂βθk +
+ ∆t 2
- − ∂2
ttmi +
- j
M i
j vβ j vγ j ∂2 βγf j eq
- + O(∆t2)
- i = 0 :
conservation of mass M 0
j ≡ 1 and the first sum is null
∂2
ttρ = −∂2 tβqβ = −∂β∂tqβ = ∂2 βγ F β γ =
- j
vβ
j vγ j ∂2 βγf j eq
and the second sum is null. Then ∂tρ + ∂β qβ = O(∆t2)
- i = α :
conservation of impulsion ∂tqα +
- j
vα
j vβ j ∂βf j eq = ∆t
- k≥3
1 sk − 1
j
vα
j vβ j (M −1)j k
- ∂βθk +
+ ∆t 2
- − ∂2
ttqα +
- j
vα
j vβ j vγ j ∂2 βγf j eq
- + O(∆t2)
ICMMES Conference, Hong Kong, July 28, 2005
From Taylor to Chapman-Enskog
−∂2
ttqα +
- j
vα
j vβ j vγ j ∂2 βγf j eq = ∂t∂βF α β +
- j
vα
j vβ j vγ j ∂2 βγf j eq
= ∂β
- j
vα
j vβ j
- ∂tf j
eq + vγ j ∂γf j eq
- = ∂β
- j
vα
j vβ j
- k
(M −1)j
k θk
= ∂β
- k≥3
j
vα
j vβ j (M −1)j k
- θk + O(∆t)
introduce the tensor Λα β
k
≡
- j
vα
j vβ j (M −1)j k
∂tqα + ∂βF α β = ∆t
- k≥3
1 sk − 1
- Λα β
k
∂βθk + ∆t 2
- k≥3
Λα β
k
∂βθk ∂tqα + ∂βF α β − ∆t
- k≥3
1 sk − 1 2
- Λα β
k
∂βθk = O(∆t2) : Chapman-Enskog !
ICMMES Conference, Hong Kong, July 28, 2005
Unidimensional heat equation with the D2Q9 LBE model
periodic boundary conditions for the y direction equilibrium velocity = 0 equilibrium energy = −2 density equilibrium (energy)2 = density relaxation of velocities : vx
∗ = vx + s1 (0 − vx)
κ = 1 3 1 s1 − 1 2
- senergy
= 1 s(energy)2 = 1 sfluxofenergy = 1, 2 sXX = sXY = 1, 4
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
∂u ∂t − ∂ ∂x
- κ∂u
∂x
- = 0 ,
0 ≤ x ≤ L , t > 0 ∂u ∂n(x = 0 , t) = Φ(t) , t > 0 ∂u ∂n(x = L , t) = 0 , t > 0 u(x , t = 0) = 0 , 0 ≤ x ≤ L unknown flux function Φ(t)
- bserved values ψ(xk, t) ≈ u(xk, t)
Inverse problem : does exists a function
- ψ(xk, •)
- k −
→ Φ(•) ? “cost” error functional : J(Φ) = 1 2
- kobserved
T
- t=0
| u(xk, t) − ψ(xk, t) |2
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
use the linearity of the problem ! Elementary discrete solution : θ(j ∆x , n ∆t) ≡ θj(n ∆t) ∂θ ∂t − ∂ ∂x
- κ ∂θ
∂x
- = 0 ,
0 ≤ x ≤ L , t > 0 ∂θ ∂n(x = 0 , t = 1) = 1 , 0 < t < ∆t ∂θ ∂n(x = 0 , t) = 0 , t > ∆t ∂θ ∂n(x = L , t) = 0 , t > 0 θ(x , t = 0) = 0 , 0 ≤ x ≤ L then Φ(t) =
N
- n=1
ϕn δ(t − n ∆t) and u(j ∆x , t) =
N
- n=1
ϕn θj(t − n ∆t) .
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
Data for Bouzidi experiment ∆x = 1 ∆t = 1 0 ≤ x ≤ L = 110 0 ≤ t ≤ T = 500 Φ(t) ≡ 0 if t ≥ N = 70
- bservation points : k = 10, 20, 30, 40, 50
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
0.1 0.2 0.3 0.4 0.5 0.6 0.7 20 40 60 80 100 grid number initial time time=50 time=100 time=150 time=200 time=250 time=300 time=350 time=400 time=450
Elementary temperature θj(n ∆t) for a unity flux
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
0.02 0.04 0.06 0.08 0.1 2 4 6 8 10 12 14 grid number (zoom near the left boundary) initial time time=50 time=100 time=150 time=200 time=250 time=300 time=350 time=400 time=450
Elementary temperature θj(n ∆t) for a unity flux
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 2 4 6 8 10 12 grid number (zoom near the left boundary) initial time time=2 time=3 time=4 time=5 time=6 time=7 time=8 time=9 time=10
Elementary temperature θj(n ∆t) for a unity flux ; first times
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
0.02 0.04 0.06 0.08 0.1 4 6 8 10 12 grid number (zoom near the left boundary) initial time time=2 time=3 time=4 time=5 time=6 time=7 time=8 time=9 time=10
the elementary temperature is always positive
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
- 0.002
0.002 0.004 0.006 0.008 0.01 0.012 0.014 4 6 8 10 12 14 grid number (zoom near the left boundary) initial time time=2 time=3 time=4 time=5 time=6 time=7 time=8 time=9 time=10
even after a big zoom !
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
5 10 15 20 25 30 35 40 20 40 60 80 100 grid number initial time time=50 time=100 time=150 time=200 time=250 time=300 time=350 time=400 time=450
“measured” temperature for a complex flux Φ(t)
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
5 10 15 20 1 2 3 4 5 6 7 8 grid number (zoom near the left boundary) initial time time=2 time=3 time=4 time=5 time=6 time=7 time=8 time=9 time=10
“measured” temperature for a complex flux Φ(t)
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
u(j ∆x , t) =
N
- n=1
ϕn θj(t − n ∆t) J(Φ) = 1 2
- k
T
- t=0
| u(xk, t) − ψ(xk, t) |2 J(Φ) = 1 2
- k
T
- t=0
|
N
- n=1
ϕn θj(t − n ∆t) − ψ(xk, t) |2 minimize the error functional J(•) relatively to the function Φ(•) δJ =
- k
T
- t=0
- N
- n=1
ϕn θj(t − n ∆t) − ψ(xk, t)
N
- m=1
δϕm θj(t − m ∆t)
- ∂J
∂ϕm =
N
- n=1
k T
- t=0
θj(t−n ∆t) θj(t−m ∆t)
- ϕn −
- k,t
ψ(xk, t) θj(t−m ∆t)
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
matrix Am n =
- k
T
- t=0
θj(t − n ∆t) θj(t − m ∆t) right hand side bn =
- k
T
- t=0
ψ(xk, t) θj(t − m ∆t) solve the linear system of order N = 70 : A • ϕ = b . What sensibility relatively to the errors of measure ?
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
5 10 15 20 100 200 300 400 500 time heat flux at the left boundary
Computed heat flux at the left boundary
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
5 10 15 20 10 20 30 40 50 60 70 time heat flux at the left boundary
Computed heat flux at the left boundary
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
Lallemand experiment : periodic boundary conditions ∂u ∂t − ∂ ∂x
- κ∂u
∂x
- = 0 ,
0 ≤ x ≤ L , t > 0 u(x = 0 , t) = u(x = L , t) t > 0 u(x = 0 , t) = 1 , 0 < t < ∆t u(x = 0 , t) = 0 , t > ∆t Two methods : LBE / Heun finite difference method ∂2u ∂x2 = five points fourth order accurate scheme Heun scheme for the temporal evolution :
- u = u(n ∆t) + ∆t ∂u
∂t
n
u((n + 1) ∆t) = u(n ∆t) + 1 2 ∂u ∂t
n
+
- ∂u
∂t
- ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
0.01 0.02 0.03 0.04 0.05 50 100 150 200 250 300 350 400 abscissa time=10 time=50 time=100 time=500
Evolution of the “elementary” temperature
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
0.01 0.02 0.03 0.04 0.05 10 20 30 40 50 abscissa (zoom near the left boundary) time=10 time=50 time=100 time=500 Gaussian reference
Evolution of the “elementary” temperature
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
correlation matrix for observation at the position k ∆x
- A(k)
- m n =
T
- t=0
θj(t − n ∆t) θj(t − m ∆t)
- rthogonal decomposition :
A = U • S • V t U and V orthogonal matrices S diagonal matrices condition number = min
j
Sj max
j
Sj
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
- 18
- 16
- 14
- 12
- 10
- 8
- 6
- 4
- 2
5 10 15 20 abscissa of the observation point kappa = 1 kappa = 2 kappa = 3 kappa = 4 kappa = 5
Logarithm of the condition number ; Boltzmann scheme
ICMMES Conference, Hong Kong, July 28, 2005
Inverse methodology for a linear thermal problem
- 18
- 16
- 14
- 12
- 10
- 8
- 6
- 4
- 2
5 10 15 20 abscissa of the observation point kappa = 1 kappa = 2 kappa = 3 kappa = 4 kappa = 5
Logarithm of the condition number ; Heun scheme
ICMMES Conference, Hong Kong, July 28, 2005
Adjoint Lattice Boltzmann Equation
- LBE discrete dynamics
f j(x, (n + 1) ∆t) = f j
∗(x − vj ∆t , n ∆t) ≡ (A • C(f, σ))(x , n ∆t)
the observation : ψ(xk, t) is supposed to be a simple function ϕ(•) of the state f(xk, t) : ϕ(f(xk, t)) ≈ ψ(xk, t)
- introduce the cost function
J(σ) ≡ 1 2
- k
t=T −∆t
- t=0
| ϕ(f(xk, t)) − ψ(xk, t) |2
- treat the LBE dynamics as a constraint (Pontryaguine, J.L. Lions) :
introduce the associated Lagrangian L ≡ J(σ) +
- x
t=T −∆t
- t=0
- p(x, t+∆t) , f(x, t+∆t) − (A • C(f, σ))(x , t)
- ICMMES Conference, Hong Kong, July 28, 2005
Adjoint Lattice Boltzmann Equation
δL = δ 1 2
- k
t=T −∆t
- t=0
| ϕ(f(xk, t)) − ψ(xk, t) |2
- +
+
- x
t=T −∆t
- t=0
δ
- p(x, t + ∆t) , f(x, t + ∆t) − (A • C(f, σ))(x , t)
=
- x
k t=T −∆t
- t=0
δx = xk
- ϕ(f(xk, t)) − ψ(xk, t)
∂ϕ ∂f δf(x, t) +
t=T −∆t
- t=0
- δp(x, t + ∆t) , f(x, t + ∆t) − (A • C(f, σ))(x , t)
- +
t=T −∆t
- t=0
- p(x, t + ∆t) , δf(x, t + ∆t) − (A • dCt(f, σ))• δf(x , t)
- −
t=T −∆t
- t=0
- p(x, t + ∆t) , A • ∂C
∂σ (f, σ)
- δσ
- ICMMES Conference, Hong Kong, July 28, 2005
Abel transform : integrate by parts the discrete sum
t=T −∆t
- t=0
- p(x, t + ∆t) , δf(x, t + ∆t) − (A • dCt(f, σ))• δf(x , t)
- =
=
t=T
- t=∆t
- p(x, t) , δf(x, t)
- −
t=T −∆t
- t=0
- p(x, t+∆t) , (A • dCt(f, σ))• δf(x , t)
- =
- p(x, T) , δf(x, T)
- +
t=T −∆t
- t=∆t
- p(x, t) , δf(x, t)
- −
t=T −∆t
- t=∆t
- dCt• At • p(x, t+∆t) , δf(x , t)
- −
- dCt• At • p(x, ∆t) , δf(x , 0)
- =
- p(x, T) , δf(x, T)
- −
- dCt• At • p(x, ∆t) , δf(x , 0)
- +
t=T −∆t
- t=∆t
- p(x, t) − dCt• At • p(x, t + ∆t) , δf(x, t)
- ICMMES Conference, Hong Kong, July 28, 2005
Variation of the Lagrangian
δL =
- x
- −
- dCt• At • p(x, ∆t) , δf(x , 0)
- +
- p(x, T) , δf(x, T)
- +
t=T −∆t
- t=∆t
- p(x, t) − dCt• At • p(x, t + ∆t) +
+
- k
δx = xk
- ϕ(f(xk, t)) − ψ(xk, t)
- ,
δf(x, t)
- +
t=T −∆t
- t=0
- δp(x, t + ∆t) , f(x, t + ∆t) − (A • C(f, σ))(x , t)
- −
t=T −∆t
- t=0
- p(x, t + ∆t) , A • ∂C
∂σ (f, σ)
- δσ
- and if the discrete dynamics LBE is satisfied,
δL ≡ δJ !!
ICMMES Conference, Hong Kong, July 28, 2005
Adjoint Lattice Boltzmann Equation
Lagrange multiplyer p(x, t) “adjoint state” integrate by parts the evaluation of δL eliminate the terms in δf in order to impose δL ≡ δJ then obtain a “backward linearized equation” for multiplyer p(x, t) also called “Adjoint Lattice Boltzmann Equation” : p(x, t) = dCt• At • p(x, t + ∆t) −
- k
δx = xk (ϕ(f(xk, t)) − ψ(xk, t)) ∂ϕ ∂f final condition p(x, T) = 0 gradient of the cost function relatively to to the parameter σ : ∂J ∂σ = −
t=T −∆t
- t=0
- p(x, t + ∆t) , A • ∂C
∂σ (f, σ)
- ICMMES Conference, Hong Kong, July 28, 2005
Gradient algorithm for optimization
minimize J(•) relatively to the parameter σ initialization : fix σ = σ0 (1) compute ∂J ∂σ :
- simulate the LBE process with a given σ and obtain f(x, t ; σ)
- simulate the ALBE process with a given f(•, • ; σ) :
final condition p(x, T) = 0 p(x, t) = dCt• At • p(x, t+∆t) −
- k
δx = xk (ϕ(f(xk, t)) − ψ(xk, t)) ∂ϕ ∂f and obtain p(x, t ; σ)
- ∂J
∂σ = −
t=T −∆t
- t=0
- p(x, t + ∆t ; σ) , A • ∂C
∂σ (f(x, t ; σ) , σ)
- evaluate a new value of σ :
σnew = σ − ρ ∂J ∂σ , ρ > 0 . go to (1)
ICMMES Conference, Hong Kong, July 28, 2005
ALBE : first results proposed by Tekitek
−14 −10 −7 −4 −1 −4.87299 −4.8729 −4.8728 −4.8727 −4.8726 −4.87256 −4.87299 −4.8729 −4.8728 −4.8727 −4.8726 −4.87256
Numerical evaluation of the cost function by finite differences
ICMMES Conference, Hong Kong, July 28, 2005
ALBE : first results proposed by Tekitek
0.1686 1 1.429 −1.9.10−6 −1.10−6 0.0 1.10−6 2.10−6 3.10−6 4.10−6 5.10−6 6.10−6 6.7.10−6 0.1686 1 1.429 −1.9.10−6 −1.10−6 1.10−6 2.10−6 3.10−6 4.10−6 5.10−6 6.10−6 6.7.10−6
Gradient of the cost function computed with finite differences (×) and with the adjoint method (+)
ICMMES Conference, Hong Kong, July 28, 2005
ALBE : first results proposed by Tekitek
10 20 30 40 −5.0 −4.0 −3.0 −2.0 −1.0 0.0 −5.0 −4.0 −3.0 −2.0 −1.0 0.0
Logarithm of the error vs iteration of the gradient method
ICMMES Conference, Hong Kong, July 28, 2005
ALBE : first results proposed by Tekitek
10 20 30 −15 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 −15 −14 −13 −12 −11 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1
Logarithm of the s5 error (top), of the s8 error (middle) and of the cost function (bottom) vs gradient iteration
ICMMES Conference, Hong Kong, July 28, 2005
ALBE : first results proposed by Tekitek
1 100 200 300 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 1 2 3 1 100 200 300 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 1 2 3