[PPT] - Stochastic Lagrangian Models Mireille Bossy Tosca Laboratory PowerPoint Presentation

SLIDE 1

Stochastic Lagrangian Models

Mireille Bossy

Tosca Laboratory Sophia Antipolis – M´ editerann´ ee

RICAM Special Semester Computational Methods in Science and Engineering School 2: Numerics for Stochastic Partial Differential Equations and their Applications December 5-9, 2016, Linz

SLIDE 2

Introduction

Stochastic Lagrangian Models (SLM) and SPDEs

In fluid mechanics, SLM are used to model the uncertainty (or the

variability) resulting from a turbulent flow. The Navier Stokes equation is an Eulerian description of the flow that gives its velocity field (U (1), U (2), U (3)) and the pressure P ∂tU (ω) + (U (ω) · ∇)U (ω) = ν△U (ω) − 1 ̺∇P(ω) + F(ω), at (t, x) t > 0, x ∈ R3, ∇ · U (t, x)(ω) = 0, t ≥ 0, x ∈ R3, U (0, x)(ω) = U0(x)(ω), x ∈ R3. (in case of an incompressible flow with constant mass density ̺)

2

SLIDE 3

Introduction

Stochastic Lagrangian Models (SLM) and SPDE

In fluid mechanics, SLM are used to model the uncertainty (or the

variability) resulting from a turbulent flow. The Navier Stokes equation is an Eulerian description of the flow that gives its velocity field (U (1), U (2), U (3)) and the pressure P ∂tU (ω) + (U (ω) · ∇)U (ω) = ν△U (ω) − 1 ̺∇P(ω) + F(ω), at (t, x) t > 0, x ∈ R3, ∇ · U (t, x)(ω) = 0, t ≥ 0, x ∈ R3, U (0, x)(ω) = U0(x)(ω), x ∈ R3. (in case of an incompressible flow with constant mass density ̺)

3

SLIDE 4

Outline Stochastic particle methods for McKean SDEs Lagrangian viewpoint in turbulence modeling

Applications Atmospheric boundary layer simulation Wind farm simulation

Well posedness results of SLM (toy models)

boundary condition mass constraint

Numerical analysis for SLM (toy models)

4

SLIDE 5

Probabilistic Interpretation of PDEs

Linear framework, forward parabolic PDE (in one dimension) A PDE of conservation form ∂ρt ∂t = 1

2σ2 ∂2ρt

∂x2 − ∂ ∂x (B(t, x)ρt) in R, ρt=0 = ρ0. (1) B(t, x), bounded measurable (1) is a Fokker-Planck Equation. Set (µ, f) =

R f(x)µ(dx).

ρ is a weak solution of (1), if for any test function f ∈ C2

b(R).

(ρt, f) = (ρ0, f) + t (ρs, 1 2σ2f ′′ + B(s, ·)f ′)ds. (2) When ρ0 is a probability measure, (1) (or (2)) describes the dynamics of the time-marginals of a stochastic process.

5

SLIDE 6

Associated martingale problem

Let (Xt)t≥0 the canonical process on (C([0, +∞), R), B(C([0, +∞), R))). For a probability P on (C([0, +∞), R), B(C([0, +∞), R))), we set Pt = P ◦ X−1

t

.

Definition - Martingale Problem

P on (C([0, +∞), R), B(C([0, +∞), R))) is a solution of the martingale problem associated to L = 1

2σ2 ∂2 ∂x2 + B(t, x) ∂ ∂x if

i) P0 = ρ0, ii) ∀f ∈ C2

b(R), the process

f(Xt) − f(X0) − t 1 2σ2f ′′(Xs) + B(s, Xs)f ′(Xs)

ds is a P-martingale.

6

SLIDE 7

Existence and uniqueness for X solution of Xt = X0 + t B(s, Xs)ds + σWt. For any test function f Ef(Xt) = Ef(X0) + t E 1 2σ2f ′′(Xs) + B(s, Xs)f ′(Xs)

ds.

(3) By setting ρt = P ◦ X−1

t

the t-time-marginal of P, Ef(Xt) =

R f(x)ρt(dx) = (ρt, f) and (3) rewrites (2) :

(ρt, f) = (ρ0, f) + t (ρs, 1 2σ2f ′′ + B(s, ·)f ′)ds.

7

SLIDE 8

Numerical approximation, particle method

Main idea : approximate the measure ρt by the empirical measure of a set of N independent trials of the r.v. Xt. Introducing N particles of independent dynamics Xi

t = Xi 0 +

t B(s, Xi

s)ds + σWi t,

i = 1, . . . , N, where ((Xi

0)i=1,...,N) are i.i.d with law ρ0, independent of the family of

independent Brownian motions ((Wi

t), i = 1, . . . , N).

U

N t = 1

N

i=1

δXi

t.

Do we have the convergence of the empirical measure U

N t toward ρt ?

8

SLIDE 9

(Basic) numerical analysis

The parabolic problem has a smooth solution

∀t > 0, ρt(dx) = ρt(x)dx. Introduce a smoothing empirical measure by a cut–off function U

N,ε t

(x) =

φε ∗ U

N t

(x) = 1

N

i=1

φε(x − Xi

t)

where φε(x) = 1

εφ( x ε) and φ : R → R smooth and such that

R φ(x)dx = 1.

At a given time t, we want the rate of convergence of the error : sup

x∈R

E

ρt(x) − U

N,ε t

(x)

.

9

SLIDE 10

Set ρε

t (x) = (ρt ∗ φε) (x).

If φ is a θ-order cut–off (

R

xqφ(x)dx = 0, ∀q ≤ θ − 1,

R

|x|θ |φ(x)| dx < ∞) and if ρt is smooth enough ρε

t (x) − ρt(x)L∞(R) ≤ Cεθ

∂θρt

∂xθ

L∞(R)

and sup

x∈R

E

ρ(t, x) − U

N,ε t

(x)

≤ Cεθ + sup

x∈R

E

ρε

t (x) − U N,ε t

(x)

with

ρε

t (x) =

R

φε(x − y)ρt(y)dy = Eφε(x − Xt) U

N,ε t

(x) = 1 N

N

i=1

φε(x − Xi

t).

The convergence is ensured by the strong law of large numbers.

10

SLIDE 11

More precisely, E

U

N,ε t

(x) − ρε

t (x)

= E
1

N

i=1
φε(x − Xi

t) − Eφε(x − Xt)

≤

1 √ N

E (φε(x − Xt) − Eφε(x − Xt))2

≤ 1 √ N

E (φε(x − Xt))2

≤ 1 √ N 1 √ε φ

1 2

L∞(R)

Eφε(x − Xt)

= 1 √ N 1 √ε φ

1 2

L∞(R)

Eρt(Yε − x)

≤ 1 √ εN φ

1 2

L∞(R) ρt

1 2

L∞(R) .

11

SLIDE 12

Adding a time discretisation of Xi

X

i (k+1)∆t = X i k∆t + ∆tB(k∆t, X i k∆t) + σ(Wi (k+1)∆t − Wi k∆t).

E

U

N,ε k∆t(x) − U N,ε,∆t k∆t

(x)

≤

1 N

N

i=1

E

φε(x − Xi

k∆t) − φε(x − X i k∆t)

≤

sup

x∈R

|φ′

ε(x)| E

X1

k∆t − X 1 t∆t

≤

C∆t ε2 (could by reduced to C∆t ε ).

Lemma

sup

x∈R

E

u(t, x) − U

N,ε,∆t t

(x)

+ E
u(t, x) − U

N,ε,∆t t

(x)

L1(R)

≤ C1εθ + C2 √ εN + C3 ∆t ε .

12

SLIDE 13

McKean Vlasov Fokker-Planck prototype

     ∂ρt ∂t = 1 2

d

i,j=1

∂2 ∂xi∂xj (aij[x, ρt]ρt) −

d

i=1

∂ ∂xi (bi[x, ρt]ρt) ρ0 a given probability measure on Rd. (4) where, for a probability measure u b[x, u] =

Rd b(x, y)u(dy); for b(x, y) Rd-valued,

a[x, u] = σ[x, u]tσ[x, u], σ[x, u] =

Rd σ(x, y)u(dy); σ(x, y) (d × k) matrix-valued.

Distribution equation : for any test function f ∈ C2

b(Rd),

(ρt, f) = (ρ0, f) + t  ρs, 1 2

d

i,j=1

aij[x, ρs] ∂2f ∂xi∂xj +

d

i=1

bi[x, ρs] ∂f ∂xi   ds (5)

13

SLIDE 14

Probabilistic point of view

Theorem [M´ el´ eard 95]

for b and σ Lipschitz on R2d, for ρ0 admitting a second order moment, there is existence and trajectorial uniqueness (given X0 ∈ L2(Ω) and (Wt)) and in law of the solution to the SDE Xt = X0 + t σ[Xs, ρs]dWs + t b[Xs, ρs]ds, where ρs is the law of Xs. (Proof from Sznitman 89, with σ = 1, and bounded b ). P2 = {probability P on C([0, T], Rd) such that EP(supt∈[0,T] |X2

t |) < ∞}, endowed

with the Wasserstein metric. Φ :P2 ∋ m → Law(Xm

t , t ∈ [0, T]) ∈ P2

Xm

t = X0 +

t σ[Xm

s , ms]dWs +

t b[Xm

s , ms]ds.

admits a fix point in P2.

14

SLIDE 15

Associated particle system

Idea : replace ρt = Law(Xt) by the empirical measure µN = 1 N

N

i=1

δXi,N. Xi,N

t

= Xi

0 +

t σ[Xi,N

s , µN s ]dWi s +

t b[Xi,N

s , µN s ]ds, i = 1, . . . , N.

(6) where ((Xi

0)i=1,...,N) are i.i.d with law ρ0, independent of the family of

independent Brownian motions ((Wi

t), i = 1, . . . , N).

(6) is well-posed when b(x, y) σ(x, y) are Lipschitz. What convergence of µN

t toward ρt ?

15

SLIDE 16

Propagation of chaos (Sznitman 89)

Let E be a separable metric space and ν a probability measure on E. A sequence of symmetric probabilities νN on EN is ν-chaotic if for any φ1, . . . , φk ∈ Cb(E; R), k ≥ 1, lim

N→∞

νN, φ1 ⊗ . . . ⊗ φk ⊗ 1 . . . ⊗ 1
=

k

l=1

ν, φl. Application : Let (X1,N, . . . , XN,N), the set of N particles satisfying (6), valued in (C([0, T]; Rd))N, with law PN.

Theorem (M´ el´ eard 95)

If the initial of (X1,N , . . . , XN,N ) is ρ0-chaotic, the chaos propagates and PN is P-chaotic where P is the unique law of the McKean SDE.

16

SLIDE 17

The tow following assertions are equivalent : (Tanaka 82)

PN is P-chaotic
µN converges weakly toward the (deterministic) measure P

Theorem (M´ el´ eard 95)

sup

N

E

sup

0≤t≤T

|Xi,N

t

|2

+ sup

N

E

sup

0≤t≤T

|Xi

t|2

< ∞

sup

N

NE

sup

0≤t≤T

|Xi,N

t

− Xi

t|2

< ∞,

where Xi

t = Xi 0 +

t σ[Xi

s, ρs]dWi s +

t b[Xi

s, ρs]ds

with ρt = P ◦ X−1

t

.

17

SLIDE 18

Numerical approximation

Xt = X0 + t σ[Xs, ρs]dWs + t b[Xs, ρs]ds, with ρs = Law(Xs). Spatial discretisation : Xi,N

t

= Xi

0 +

t σ[Xi,N

s , µN s ]dWi s +

t

)

b[Xi,N

s , µN s ]ds, i ≤ N.

Time-Discretisation : ∆t > 0, such that T = K∆t for K ∈ N,                  ¯ Xi,N

(k+1)∆t = ¯

Xi,N

k∆t + ( 1

N

j=1

σ(X

i,N k∆t, X j,N k∆t))(Wi (k+1)∆t − Wi k∆t)

+ ∆t( 1 N

N

j=1

b(¯ Xi,N

k∆t, ¯

Xj,N

k∆t)),

¯ Xi,N

t=0 = Xi 0.

U

N,ε,∆t k∆t

(x) = 1 N

N

i=1

φε(x − ¯ Xi,N

k∆t), pour un cut-off donn´

e φ.

18

SLIDE 19

Rate of convergence

Theorem (Bossy & Talay 97)

Dimension d = 1. Assume that ρ0(·) ∈ C2

K(R), b(·, ·) ∈ C2,2 b (R2) and s(·, ·) ∈ C3 b(R2).

Assume the strong ellipticity of the differential operator L. Then ρt(·) ∈ W2,1(R) for all t ∈ [0, T]. For a 2-order cut-off φ(x), ∀ k ≤ K, EU

N,ε,∆t k∆t

− ρk∆tL1(R) ≤ Cb,σ,T,ρ0

1

ε √ N + ε2 + √ ∆t ε

.

H(x) := 1{x≥0}, Heaviside function E(H ∗ ρk∆t)(x) − 1 N

N

i=1

H(x − X

i,N k∆t)L1(R) ≤ Cb,σ,T,ρ0

1 √ N + √ ∆t

.

19

SLIDE 20

Rate of convergence

Theorem (Antonelli & Kohatsu-Higa 02)

d = 1. b(·, ·) and s(·, ·) C∞(R2), with bounded derivatives. A non degenerescency condition for L (Malliavin calculus). φ(x) is the Gaussian density on R. ε2 = ∆t. EU

N,ε,∆t k∆t

− ρk∆tL1(R) ≤ Cb,σ,T,ρ0

1

√ ∆tN + ∆t + 1 √ N

.

H(x) := 1x≥0, Heaviside function, E(H ∗ ρk∆t)(x) − 1 N

N

i=1

H(x − X

i,N k∆t)L1(R) ≤ Cb,σ,T,ρ0

1 √ N + ∆t

.

20

SLIDE 21

Incompressible 2d-Navier Stokes equation

V(t, x) = (v1(t, x), v2(t, x)) and p(t, x) solution to    ∂tV = ν∆V − (V · ∇)V − ∇p, t ≥ 0, x ∈ R2, ∇ · V = ∂x1v1 + ∂x2v2 = 0, V(0, x) = V0(x). (7) Vorticity formulation : U = ∇ × V = ∂x1v2 − ∂x2v1 obtained by derivation.      ∂U ∂t (t, x) = ν∆U(t, x) − ∇.

U(

K ∗ U)(t, x)

, x ∈ R2,
V(x, t) = (

K ∗ U)(x, t), U(0, x) = (∇ × V0) (x), (8) with

K(x1, x2) =

1 2π(x2

1 + x2 2)(−x2, x1).

Chorin 73 : Random Vortex Method.

21

SLIDE 22

N. S. incompressible 2d : probabilistic Interpretation

Marchioro & Pulvirenti 82, Osada 85, M´ el´ eard 01

Introduce a martingale problem, where the weighted marginal density

laws ( Pt)t≥0 from a solution P, solve the vortex equation and V(t, x) = ( K ∗ Pt)(x).

Convergence of the empirical measure of the associated particle system

to the unique solution of the martingale problem (chaos propagation).

Convergence of the velocity field given by the empirical measure to the

solution of the NS equation VN(t, x) = ( K ∗ µN

t )(x) to V(t, x).

22

SLIDE 23

Numerical analysis : by splitting method

Convergence of the vortex algorithm, Beale & Majda 82, Goodman 87, . . .
A first rate of convergence result Long 88:

Assume that U0 has compact support. Fix h, Λh = {h.i, i ∈ Z2}. Uh

0(x) =

Λh∩Supp(U0)

h2wiδ(x−αi) φ ∈ CL(R2) cut–off function of order m, with fast decay at infinity The particle system : dXε,h,i

t

=

j

wjh2Kε

Xε,h,i

t

− Xε,h,j

t

dt +

√ 2νdWi

t.

Vε,h(t, x) =

j

wjh2Kε

x − Xε,h,j

t

.

23

SLIDE 24

N. S. incompressible 2d: rate of convergence ?

(Long 88) Choose ε ≥ hq, q < 1. Assume that V(t, x) is smooth enough. sup

0≤t≤T

Vε,h(t, x)(ω) − V(t, x)
Lp(B(R0)) ≤ C
εm +

h ε L ε + h| ln h|

except for ω in A, an event of probability smaller than hC′C,

with C > 0 and C′ > 0 depending in all the parameters (T, L, m, p, q, R0, V0, ∂i,jV).

24

SLIDE 25

Incompressible N.S. in a bounded domain O ⊂ R2

V(t, x) = (v1(t, x), v2(t, x)) and p(t, x) solve    ∂tV = ν∆V − (V · ∇)V − ∇p, t ≥ 0, x ∈ O, ∇ · V = ∂x1v1 + ∂x2v2 = 0 V(0, x) = V0(x), V(t, x) = 0, x ∈ ∂O. (9) Vorticity equation : U = ∂x1v2 − ∂x2v1      ∂U ∂t (t, x) = ν∆U(t, x) − ∇. (UV(t, x)) , x ∈ O, ∇ · V = 0 ⇒ V = ∇⊥ψ = (∂x2ψ, −∂x1ψ) with U = −∆ψ, Ut=0 = U0 = ∇ × V0, Choose ψ(t, x) = −

O γ(x, y)U(t, y)dy, where γ is the Green function of ∆ on O

with homogeneous Dirichlet condition V(t, x) = −

O

∇⊥

x γ(x, y)U(t, y)dy.

Resulting boundary condition for the vorticity : ∀x ∈ ∂O, ψ(t, x) = 0 ⇒ ∂τψ = 0 ⇒ V.n = 0.

25

SLIDE 26

Incompressible N.S. in O : probabilistic interpretation

We need to vanish the tangential velocity V.τ = 0 on ∂O, from the vorticity created at the boundary (Vortex sheet method, Chorin 73, 78)

Benachour Roynette Vallois 01 : Branching process at the boundary for

the Vorticity SDE Difficult to simulate the associated particle system.

Jourdain & M´

el´ eard 04 : give a probabilistic interpretation of the Vorticity equation with additional Neumann boundary condition      ∂U ∂t (t, x) = ν∆U(t, x) − ∇. (UV(t, x)) , x ∈ O, Ut=0 = U0 = ∇ × V0 inO, ∂nU(t, x) = g(t, x) in ]0, T[×∂O.

Cottet 89 proposes a non linear Neumann boundary condition, that

exactly produces the no slip boundary condition of N.S.

26

SLIDE 27

Outline Stochastic particle methods for McKean SDEs Lagrangian viewpoint in turbulence modeling

Applications Atmospheric boundary layer simulation Wind farm simulation

Well posedness results of SLM (toy models)

boundary condition mass constraint

Numerical analysis for SLM (toy models)

27

SLIDE 28

SLM for Atmospheric Boundary Layer simulation

28

SLIDE 29

Multiphase flows

29

SLIDE 30

Statistical approach of turbulent flows

Eulerian description

◮ The incompressible Navier Stokes equation in R3, for the velocity field

(U (1), U (2), U (3)) and the pressure P, with constant mass density ̺ ∂tU (ω) + (U (ω) · ∇)U (ω) = ν△U (ω) − 1 ̺∇P(ω), t > 0, x ∈ R3, ∇ · U (ω) = 0, t ≥ 0, x ∈ R3, U (0, x)(ω) = U0(x)(ω), x ∈ R3.

◮ Ensemble averages as expectations:

U (t, x) :=

Ω

U (t, x, ω)dP(ω). Reynolds decomposition U (t, x, ω) = U (t, x) + u′(t, x, ω), P(t, x, ω) = P(t, x) + p′(t, x, ω) The random field u′(t, x, ω) is the turbulent part of the velocity.

30

SLIDE 31

Reynolds-Averaged Navier-Stokes (RANS) Equation

Assuming Reynolds decomposition, with constant mass density ̺, we get ∂tU (i) +

3

j=1

U (j)∂xjU (i) +

3

j=1

∂xju′

iu′ j = ν△U (i) − 1

̺∂xiP, ∇ · U = 0, t ≥ 0, x ∈ R3, U (0, x) = U0(x), x ∈ R3, closed by a chosen parameterization or a model equation of the Reynolds stress tensor (u′

iu′ j = U (i)U (j) − U (i)U (j), i, j).

Turbulent kinetic energy tke(t, x) := 1 2

3

i=1

u′

iu′ i(t, x)

Pseudo-dissipation ε(t, x) := ν

3

i=1

3

k=1

∂xku′

i∂xku′ i(t, x)

31

SLIDE 32

Equation for the Reynolds stress (u′

iu′ j, i, j) ∂tu′

iu′ j +

U · ∇xu′

iu′ j

+

3

k=1

∂xku′

iu′ ju′ k

= −1 ̺u′

j∂xip′ + u′ i∂xjp′

velocity pressure gradient tensor Πij

+ν△xu′

iu′ j

+ ν

3

k=1

∂xku′

i∂xku′ j

dissipation tensor εij

−

3

k=1
u′

iu′ k∂xkU (i) + u′ ju′ k∂xkU (j)

turbulence production tensor Pij

32

SLIDE 33

Alternative viewpoint

The PDF approach

Let ρEuler(t, x; V) be the probability density function of the (Eulerian) random field U (t, x), then U (i)(t, x) =

R3 V(i)ρEuler(t, x; V)dV,

U (i)U (j)(t, x) =

R3 V(i)V(j)ρEuler(t, x; V)dV.

The closure problem is reported on the PDE satisfied by the probability density function ρEuler. This is the so-called PDF method. In his seminal work, Stephen B. Pope proposes to model the PDF ρEuler with the Lagrangian probability density function,

r equivalently with a Lagrangian description of the flow.

33

SLIDE 34

Lagrangian description for turbulent flow

On a given probability space (Ω, F, P), consider the fluid particle state vector (Xt, Ut, θt) satisfying dXt =Utdt, dUt =

−1

̺∇xP(t, Xt) − G(t, Xt) (Ut − U (t, Xt)) + Ft

dt

+

C(t, Xt)ε(t, Xt)dWt,

dθt =D1(t, Xt, θt)dt + D2(t, Xt, θt)d Wt. (W, W) is a (3d × 1d)-Brownian motion. U (i)(t, x), U (i)U (j)(t, x) are computed as conditional expectations. ε, C, G, D1, D2 are determined by a chosen Reynolds-averaged Navier-Stokes closure.

34

SLIDE 35

Compute the Reynolds averages U (i) and U (i)U (j)

Let ρLagrangian(t; x, V, ψ) be the probability density function of (Xt, Ut, θt). Eulerian PDF versus Lagrangian PDF:

(case of incompressible flow, constant mass of particles)

ρEuler(t, x; V, φ)dVdφ = ρLagrangian(t; x, V, ψ)

R4 ρLagrangian(t; x, U, ψ)dUdψ dVdφ

For any random variable F(Ut, θt), with finite moment, F(U , ϑ)(t, x) =

R4 F(V, ψ)ρEuler(t, x; V, ψ)dVdψ

= E

F(Ut, θt)
Xt = x
In particular,

U (i)(t, x) = E

U(i)

t

Xt = x
U (i) U (j)(t, x) = E
U(i)

t

U(j)

t

Xt = x
35

SLIDE 36

Fluid particle model family

On a given probability space (Ω, F, P), consider the fluid particle state vector (Xt, Ut, ψt) satisfying dXt =Utdt dUt =

−1

̺∇xP(t, Xt) − G(t, Xt) (Ut − U (t, Xt)) + Ft

dt

+

C(t, Xt)ε(t, Xt)dWt,

dθt =D1(t, Xt, θt)dt + D2(t, Xt, θt)d Wt. (W, W) is a (3d × 1d)-Brownian motion. The Eulerian fields U (i)(t, x), U (i)U (j)(t, x), ϑU (j)(t, x), and their derivatives, are nonlinear McKean terms in this SDE. One needs to specify ε, C, G, D1, D2 imposed by the chosen turbulence closure compute ∇xP introduce boundary conditions

36

SLIDE 37

PDF description of flow via the associated Fokker Planck Equation

P ((Xt, Ut, ψt) ∈ dxdvdφ) := ρLagr.(t; x, v, φ)dxdvdφ. ∂tρLagr. + (v · ∇xρLagr.) = 1 ̺ (∇xP(t, x) · ∇vρLagr.) − ∇v · (G(t, x)(U (t, x) − v)ρLagr.) + 1 2C(t, x)ε(t, x)△vρLagr. − ∂φ(D1(t, x, φ)ρLagr.) + 1 2∂2

φφ (D2(t, x, φ)ρLagr.)

Integrating w.r.t. dvdφ gives the equation for the conservation of mass, ̺(t, x) =

ρLagr.(t; x, v, φ)dvdφ

∂t

ρLagr.dvdφ + ∇x ·
vρLagr.dvdφ
ρLagr.dvdφ
ρLagr.dvdφ
= 0 ⇐

⇒ ∂t̺ + ∇x · (̺U ) = 0. Integrating w.r.t. vdvdφ gives the RANS equation. Integrating w.r.t. vvtdvdφ gives the Reynolds stress closure, according to C, G, D1, D2.

37

SLIDE 38

The Reynolds tensor models in SLM

Simple Langevin Model see e.g [Pope 95, Minier Peirano 01, Pope 03] Durbin Speziale 94, Pope 94, Dreeben Pope 98,Durbin 93, Waclawczyk Pozorski Minier 04,      dXt = Utdt, with Ut = (u(i)

t , i = 1, 2, 3) and u′ i(t) = u(i) t

− u(i)

t

du(i)

t

= −∂xiP(t, Xt)dt +

j

Gij

u(j)

t

− u(j)

(t, Xt)dt +
C0ε(t, Xt)dB(i)

t

Gij = −CR 2 ε tkeδij +C2 ∂u(i) ∂xj C0ε = 2 3 (CRε+C2P − ε) , P = 1

2 (P11 + P22 + P33), the turbulent production term is computed with

Pij := −

k
u′

iu′ k∂u(i)

∂xk + u′

ju′ k∂u(j)

∂xk

.

ε is closed with mixing length or turbulent viscosity parametrizations.

38

SLIDE 39

The Reynolds tensor models in SLM

Corresponding Reynolds Stress Equation :

∂tu′

iu′ j +

U · ∇xu′

iu′ j

+

3

k=1

∂xku′

iu′ ju′ k =

−

3

k=1
u′

iu′ k∂xkU (i) + u′ ju′ k∂xkU (j)

− 2

1 2 + 3 4 C0 ε tke u′

iu′ j + C0εδij

versus Exact Reynolds Stress Equation :

∂tu′

iu′ j +

U · ∇xu′

iu′ j

+

3

k=1

∂xku′

iu′ ju′ k =

−

3

k=1
u′

iu′ k∂xkU (i) + u′ ju′ k∂xkU (j)

− 1

̺ u′

j∂xip′ + u′ i∂xjp′ + ν 3

k=1

∂xku′

i∂xku′ j

The dissipation tensor is reduced to the relation ν

3

l=1

∂xlu′

i∂xlu′ j = 2

3εEδij The turbulent pressure terms are closed according to the Rotta’s second order closure model : 1 ̺

p′∂xju′

i + p′∂xiu′ j

= −
1 + 3

2C0 ε tkeu′

iu′ j + 2

3CRεδij, for

1 + 3

2C0

denoting the prescribed Rotta’s constant.

39

SLIDE 40

Reynolds tensor models in SLM

IP model for anisotropic turbulence see e.g [Pope 95, Minier Peirano 01, Pope 03, Durbin Speziale 94, Pope 94, Dreeben Pope 98]Durbin 93, Waclawczyk Pozorski Minier 04,      dXt = Utdt, with Ut = (u(i)

t , i = 1, 2, 3) and u′ i(t) = u(i) t

− u(i)

t

du(i)

t

= −∂xiP(t, Xt)dt +

j

Gij

u(j)

t

− u(j)

(t, Xt)dt +
C0ε(t, Xt)dB(i)

t

Gij = −CR 2 ε tkeδij +C2 ∂u(i) ∂xj C0ε = 2 3 (CRε+C2P − ε) , P = 1

2 (P11 + P22 + P33), the turbulent production term is computed with

Pij := −

k
u′

iu′ k∂u(i)

∂xk + u′

ju′ k∂u(j)

∂xk

.

ε is closed with mixing length or turbulent viscosity parametrizations.

40

SLIDE 41

Reynolds tensor models in SLM

see e.g [Pope 95, Minier Peirano 01, Pope 03, Durbin Speziale 94, Pope 94, Dreeben Pope 98, Durbin 93, Waclawczyk Pozorski Minier 04]

       dXt = Utdt, with Ut = (u(i)

t , i = 1, 2, 3) and u′ i(t) = u(i) t

− u(i)

t

du(i)

t

= −∂xiP(t, Xt)dt +  

j

Gij

u(j)

t

− u(j)



 (t, Xt)dt +

C0ε(t, Xt)dB(i)

t

Gij = −γij − 1 2 ε tkeδij C0ε = 2 3(γij)u′

iu′ j

γij = (1 − αtke)γ

wall

ij + αtkeγ

homogeneous

ij

Elliptic blending coefficent α solves: L2∇2α − α = − 1 tke (L length scale), and −γ

homogeneous

ij

= −1 2(CR − 1) ε tkeδij + C2 ∂u(i) ∂xj −γ

wall

ij

= −7.5 ε tkeninj where n is the wall-normal unit vector.

41

SLIDE 42

Dissipation closure in the atmospheric boundary layer

collaboration with P . Drobinski (IPSL)

With notation Ut = (ut, vt, wt) and also Ut − U(t, Xt) = (u′

t, v′ t, w′ t)

Closure relation classically used in meteorology (see e.g.

Cuxart et al. 2000):

ε(t, x, y, z) = Cε tke3/2(t, x, y, z) ℓm(z) . Relation based on turbulent viscosity parametrization (O’Brien 70, Heinz 97) νturb =

u′w′

2+v′w′ 2

∂zu

2+

∂zv

2

ε = tke2 3 1 1

2 + 3 4C0

1

νparam

turb

ℓm(z) = κz1[0,zc] + κzc1[zc,750]

42

SLIDE 43

Wall-boundary condition

The velocity is reflected at level zmirror (Drebeen Pope 98, Minier

Pozorski 99)

win = −wout uin = uout − 2u′w′ w′2 wout, vin = vout − 2v′w′ w′2 wout.

Cell center

Cells {z = zc} Outward particle

Uout

Uin

Mirror particle
Inward normal

{z = zmirror} {z = z0} {z = 0} Ground

We compute u′w′ and v′w′ using the MESO-NH-style formula u′w′(t, x, y) = − u2

∗u

u2 + v2 (t, x, y),

v′w′(t, x, y) = − u2

∗v

u2 + v2 (t, x, y)

with u∗(t, xc, yc) = κ

u2(t, xc, yc, zc) + v2(t, xc, yc, zc)

log (zc/z0)

(see e.g. [Schmid and Schumann 89, Lafore et al 98, Carlotti 02])

43

SLIDE 44

Validation in neutral atmosphere

collaboration with P . Drobinski (IPSL)

Comparison with a LES method [MESO-NH use case in Drobinski et al. 2007]. Computational domain : 3000 m × 1000 m × 750 m Renormalized variances and covariances u′u′, v′v′, w′w′, u′w′, v′w′, u′v′,

44

SLIDE 45

Numerical method : stochastic particle algorithm

The PIC method

The computational space is divided in cells We use a Particles in cell (PIC) technique to compute the Eulerian fields (mean fields) U (i)(t, x), U (i)U (j)(t, x) at the center of each cell : We introduce Np particles (Xk,Np

t

, Uk,Np

t

) in D. Each cell C contains Npc particles by constant mass density constraint. K(., x) = 1(., C(x)). F(U )(t, x) ≃ 1 Np

Np

k=1

F

Uk,Np

t

K(Xk,Np

t

, x) 1 Np

Np

k=1

K(Xk,Np

t

, x)

Np

k=1

K(Xk,Np

t

, x) = Npc

45

SLIDE 46

The interacting particles system

For j = 1, . . . , Np            dXj,Np

t

= Uj,Np

t

dt, dU(i),j,Np

t

= −1 ̺ ∂P ∂xi (t, Xj,Np

t

)dt +DU(t, Xj,Np

t

)dt + BU(t, Xj,Np

t

)dW(i),j,Np

t

+jump terms for boundary condition, i ∈ {1, 2, 3}.

The coefficients DU, BU depend on the particles approximations of U ,

U (i)U (j), and their derivatives.

−1

̺ ∂P ∂xi (t, Xj,Np

t

) ensures that ∇ · U = 0 and maintains the mass density constant. ∇2P = −∂2U (i)U (j) ∂xi∂xj

46

SLIDE 47

Algorithm

Fractional step method : n∆t − → (n + 1)∆t (Pope 85, Minier Pozorski 99) Step 1 : For n∆t ≤ t ≤ (n + 1)∆t, (Xj,Np

n∆t, U(i),j,Np n∆t

, j = 1, . . . , Np), given          d Xj,Np

t

= Uj,Np

t

dt, d U(i),j,Np

t

= − 1

̺ ∂P ∂xi

(t, Xj,Np

t

)dt +D

U(t, Xj,Np t

)dt + B

U(t, Xj,Np t

)dW(i),j,Np

t

+jump terms, i ∈ {1, 2, 3} Step 2: Correction of the particles positions Xj,Np

(n+1)∆t −

→ Xj,Np

(n+1)∆t

that maintains the (discrete) uniform distribution, by solving a (discrete) optimal transport problem. Step 3: Correction of the particles velocity Uj,Np

(n+1)∆t −

→ Uj,Np

(n+1)∆t

such that ∇ · U (n+1)∆t = 0, as in Chorin-Temam projection method :

△p =

1 ∆t∇ ·

U (n+1)∆t, x ∈ D,

∂p ∂nD = 0, on ∂D

Uj,Np

(n+1)∆t =

Uj,Np

(n+1)∆t − ∆t∇p(Xj,Np (n+1)∆t) and U (n+1)∆t =

U (n+1)∆t − ∆t∇p.

47

SLIDE 48

Turbine wake effects

(see e.g. Hansen 08, R´ ethor´ e, Sørensen, Bechmann 09, Port´ e-Agel et al 10)

48

SLIDE 49

Actuator disc methods in the Lagrangian setting

dXt =Utdt dUt =

−1

ρ∇xP(t, Xt)

dt

− G(t, Xt)

Ut − U(t, Xt)
dt + C(t, Xt)dWt

+ f (t, Xt, Ut) dt + fnacelle (t, Xt, Ut) dt + fmast (t, Xt, Ut) dt. f(t, Xt, Ut) the body forces that the blades exert on the flow. Supplementary terms fnacelle(t, Xt, Ut) and fmast(t, Xt, Ut) represent the impact of the mill nacelle and mast. (B., Espina, Morice, Paris, Rousseau 14)

49

SLIDE 50

Forces

Non rotating, uniformly loaded actuator disc model: knowing the axial induction factor a fx(t, Xt, Ut) = − 1 ∆x 2a 1 − a

E [ut|Xt ∈ C]
2ex.

Rotating actuator disc :

fx(t, Xt, Ut) = −1{Xt∈C} Nb 4πr∆xUrelat(Xt, Ut)2c(r(Xt)) (CL(α) cos(φ) + CD(α) sin(φ)) (Xt, Ut), fθ(t, Xt, Ut) = 1{Xt∈C} Nb 4πr∆xUrelat(Xt, Ut)2c(r(Xt)) (CL(α) sin(φ) − CD(α) cos(φ)) (Xt, Ut).

50

SLIDE 51

Validation of with wind tunnel measures

Figure 1: Comparison of vertical profiles of u at different downstream positions x from the turbine (-1, 2, 5, 7 and 10 diameters respectively). The profiles are centred to the hub y position. The blue curve represents the wind tunnel measures, the red curve represents SDM simulation with the Rotating Actuator Disc mill model.

51

SLIDE 52

Validation with wind tunnel measures

Figure 2: Turbulence intensity with SDM, using the inflow mean velocity Uhub at the hub height given in Wu-Porte Agel 2011: I =

2/3 tke/Uhub

tke isosurface, for a (3 blades) actuator line model.

52

SLIDE 53

Figure 3: Domain for the wind tunnel scale simulations and numerical parameters

53

SLIDE 54

Wind farm evaluation

54

SLIDE 55

Wind farm evaluation

55

SLIDE 56

Outline Stochastic particle methods for McKean SDEs Lagrangian viewpoint in turbulence modeling

Applications Atmospheric boundary layer simulation Wind farm simulation

Well posedness results of SLM (toy models)

boundary condition mass constraint

Numerical analysis for SLM (toy models)

56

SLIDE 57

Well posedness of a simplified Langevin model

Consider    Xt = X0 + t

0 Usds,

Ut = U0 + t

0 B [Xs, Us; ρs] ds +

t

0 σ(s, Xs, Us) dWs,

P((Xt, Ut) ∈ dxdu) = ρt(x, u)dxdu, t ∈ [0, T], where B[x, u; γ] =   

Rd b(u, v)γ(x, v)dv
Rd γ(x, v)dv

, if

Rd γ(x, v)dv = 0,

0, elsewhere If b : R2d → Rd is bounded, by the Girsanov theorem, any weak solution (Xt, Ut, t ∈ [0, T]) has a strictly positive density (ρt(x, u), t ∈ [0, T]), and E [b(u, Ut)|Xt = x] = B[x, u; ρt].

57

SLIDE 58

Well posedness of a simplified Langevin model

◮ σ is bounded and strongly elliptic: for a := σσ∗, there exists λ > 0 s.t. for

all t ∈ (0, T], x, u, v ∈ Rd, |v|2 λ ≤

d

i,j=1

a(i,j)(t, x, u)vivj ≤ λ |v|2 .

◮ For all 1 ≤ i, j ≤ d, a(i,j) is H¨

lder continuous in the following sense:

α ∈ (0, 1] |a(i,j)(t, x, u) − a(i,j)(s, y, v)| ≤ C(|t − s|

α 2 + |x − y − v(t − s)|α + |u − v|α).

Theorem Bossy, Jabir, Talay 11

Let b ∈ Cb(R2d, Rd), let (X0, U0) s.t. EP

X0Rd + U02

Rd

< +∞. On the

previous hypotheses on the velocity diffusion coefficient σ, the system has a unique weak solution.

58

SLIDE 59

The smoothed system in the space variables

Xε

t = X0 +

t

0 Uε s ds,

Uε

t = U0 +

t

0 Bε [Xε s , Uε s ; ρε s ]ds +

t

0 σ(s, Xs, Us) dWs,

where Law(Xε

t , Uε t ) = ρε t (x, u)dxdu,

and ∀ non–negative γ in L1(R2d), Bε [x, u; γ] =

R2d b(v, u)φε(x − y)γ(y, v) dy dv
R2d φε(x − y)γ(y, v) dy dv + ε ,

for a given regularization φε ∈ C1

c(Rd) of the Dirac mass in Rd.

             dXε

t = Uε t dt,

dUε

t =

E [b(v, Uε

t )φε(x − Xε t )]

x=Xε

t ,v=Uε t

E [φε(x − Xε

t )]

x=Xε

t

+ε dt + t σ(s, Xε

s , Uε s ) dWs,

t ∈ [0, T]. (10)

59

SLIDE 60

Existence of the smoothed system

On (Ω, F, (Ft; t ∈ [0, T]) , Q),

Wi

t; t ∈ [0, T] ; N

independent Brownian motions in Rd,

{(Xi

0, Ui 0); i ∈ N} i.i.d, independent of the Brownian family, such that

Q

(X1

0, U1 0) ∈ dx du

= ρ0(x, u) dx du.

                 Xi,ε,N

t

= Xi

0 +

t

0 Ui,ε,N s

ds, Ui,ε,N

t

= Ui

0 +

t

N

j=1, j=i

b(Ui,ε,N

s

, Uj,ε,N

s

)φε(Xi,ε,N

s

− Xj,ε,N

s

) N

j=1, j=i

φε(Xi,ε,N

s

− Xj,ε,N

s

) + ε

ds + Wi

t,

i = 1, . . . , N. The sequence {πN = Law

1

N

i=1 δ{Xi,ε,N,Ui,ε,N,Wi}

, N ∈ N} is tight on

P(C([0, T]; R3d)). The sequence {Law(Xi,ε,N, Ui,ε,N, 1 ≤ i ≤ N), N} is Law(Xε, Uε)–chaotic.

60

SLIDE 61

Convergence ε → 0, uniqueness results

The linear Langevin system (I)

For (y, v) ∈ R2d, consider (Ys,y,v

t

, Vs,v

t ; t ≥ s ≥ 0) solution of the Langevin

equation        Ys,y,v

t

= y + t

s

Vs,y,v

θ

dθ, Vs,y,v

t

= v + t

s

σ(θ, Ys,y,v, Vs,y,v

θ

) dWθ. (11)

Proposition

There exists a unique weak solution to (11). In addition this solution admits a density Γ(s, y, v; t, x, u) w.r.t. Lebesgue measure such that: (i) For all (t, x, u) ∈ (0, T] × R2d, 1 ≤ i, j ≤ d, the derivatives ∂viΓ(s, y, v; t, x, u) exist and are continuous for all (s, y, v) ∈ R × R2d such that (s, y, v) = (t, x, u).

(Ultraparabolic operators of Kolmogorov type (Di Francesco&Pascucci 05))

61

SLIDE 62

Convergence ε → 0, uniqueness results

The linear Langevin system (II) Proposition (following)

(ii) Let f ∈ Cb(Rd). The function Gt,f defined by Gt,f (s, y, v) = E (f(Ys,y,v

t

, Vs,y,v

t

)) , is the unique classical solution of the Cauchy problem:    ∂sGt,f + A (s, Gt,f ) = 0 on [0, t) × R2d, lim

s→t− Gt,f (s, y, v) = f(y, v) ∀(y, v) ∈ R2d,

A(ψ)(s, y, v) = (u · ∇xψ(s, y, v)) + 1 2

d

i,j=1

(σσ∗)(i,j)(s, x, v)∂2

vi,vjψ(s, y, v)

(iii) There exists C > 0 depending only on T and λ such that sup

(y,v)∈R2d

Rd |∇vΓ(s, y, v; t, x, u)| dx du ≤

C √t − s, ∀ 0 ≤ s < t ≤ T.

62

SLIDE 63

Convergence ε → 0, uniqueness results

(S∗

t,s) is the adjoint of the transition semigroup of (Ys,y,v t

, Vs,y,v

t

), Define S′

t,s : L1(R2d; Rd) → L1(R2d; R) by

S′

t,s(h)(x, u) =

R2d (∇vΓ(s, y, v; t, x, u) · h(y, v)) dy dv.

Any solution of the (smoothed) Lagrangian model has a positive marginal density satisfying the mild equation in L1(R2d) ρt = S∗

t,0(µ0) +

t S′

t,s (ρs(·)B [· ; ρs]) ds

ρǫ

t

= S∗

t,0(µ0) +

t S′

t,s (ρǫ s(·)Bǫ [· ; ρǫ s]) ds

For all t ∈ (0, T], h, δ ∈ Rd, it holds that lim

|h|,|δ|→0 lim sup ǫ→0+

R2d |ρǫ

t (x + h, u + δ) − ρǫ t (x, u)| dx du = 0.

63

SLIDE 64

Definition - : The non linear Martingale Problem (MPε) :

A probability Pε on C([0, T]; R2d) is a weak solution of (10), or equivalently, a solution to the martingale problem (MPε) if (i) Pε ◦ (x0, u0)−1 = µ0. (ii) For all t ∈ (0, T], the time marginal Pε ◦ (xt, ut)−1 has a density ρǫ

t w.r.t.

Lebesgue measure on R2d. (iii) For all f ∈ C2

b(R2d), the process

f(xt, ut) − f(x0, u0) − t Aǫ

ρǫ

s f(s, xs, us) ds

(12) is a Pε–martingale

Any limit point π∞ of {πN, N ∈ N} has a full measure on the set of solutions
f (MPε).
Consequently: Any subsequence of {πN, N ∈ N} converges to δPε or

equivalently Law

Xi,ε,N, Ui,ε,N, Wi, i = 1, . . . , N
is Pε-chaotic.

64

SLIDE 65

The limit ε → 0

Lemma : The limit Martingale Problem (MP)

There exists P on C([0, T]; R2d) such that P is a solution to the martingale problem (MP) : (i) P ◦ (x0, u0)−1 = µ0. (ii) For all t ∈ (0, T], the time marginal P ◦ (xt, ut)−1 has a positive density ρt w.r.t. Lebesgue measure on R2d. (iii) For all f ∈ C2

b(R2d), the process

f(xt, ut) − f(x0, u0) − t Aρsf(s, xs, us) ds is a (MP)–martingale. Main difficulty : for any bounded Fs-measurable r.v Φ, lim

ε→0+ EPε

Φ

t

s

(Bε [xθ, uθ; ρε

θ] · ∇uf(xθ, uθ)) dθ

= EP
Φ

t

s

(B [xθ, uθ; ρθ] · ∇uf(xθ, uθ)) dθ

.

65

SLIDE 66

Compute the wind at a finer scale

Downscaling method for local wind assessment, based on a stochastic Lagrangian approach

66

SLIDE 67

The downscaling problem

0◦ 4◦E 8◦E 41◦N 44◦N

Mer M´ editerran´ ee France

Comparison WRF and SDM

n flat terrain;

horizontal visualization of the u component of the wind at the altitude 400 m ; time is T0 + 4 days.

67

SLIDE 68

Numerical framework

Our computational domain D corresponds to a given set of cell of a NWP solver. Top and lateral (in red) : mean-Dirichlet boundary condition ∀x ∈ ∂D, U (t, x) = Vext(t, x) where Vext is the WRF forcing. Bottom boundary condition (in blue) necessities to introduce a wall law model.

68

SLIDE 69

The boundary condition in the PDF approach

∀x ∈ ∂D, U (t, x) = Vext(t, x).

Rd vρℓ(t, x, v)dv
Rd ρℓ(t, x, v)dv = Vext(t, x).
Rd vρℓ(t, x, v, )dv =
Rd Vext(t, x)ρℓ(t, x, v)dv

⇔

Rd vρℓ(t, x, v, )dv =
Rd vρℓ(t, x, v + 2(Vext(t, x) − v))dv

⇑ ρℓ(t, x, v, ) = ρℓ(t, x, v + 2(Vext(t, x) − v)), ∀v ∈ Rd leads to specular boundary condition with jump on ∂D for ρℓ ...

69

SLIDE 70

Stochastic Downscaling Method (SDM)

Let D be an open set of R3, and a velocity Vext given at ∂D: Xt =X0 + t Usds, dUt =U0 + t

−1

̺∇P(s, Xs) − 1 2 + 3 4C0 ε(s, Xs) tke(s, Xs) (Us − U (s, Xs))

ds

+ t

C0ε(s, Xs)dWs+
0<s≤t

2 (Vext(s, Xs) − Us−) 1{Xs∈∂D}. The jump term should ensure that U (t, x) = Vext(t, x), ∀x ∈ ∂D. Generic nonlinear SDE : Xt = X0 + t

0 Us ds,

Ut = U0 + t

0 B[Xs; ρ(s)]ds +

t

0 Σ[Xs; ρ(s)]dWs−

0<s≤t

2 (Vext − Us−) 1{Xs∈∂D}, ρ(t) is the probability density of (Xt, Ut) for all t ∈ (0, T], B[x; f] =

Rd b(v)f(x, u)du
Rd f(x, u)du

1{

Rd f(x,u)du=0}, Σ[x; f] =
Rd σ(u)f(x, u)du
Rd f(x, u)du

1{

Rd f(x,u)du=0}.

70

SLIDE 71

Well-posedness of the SLM with the no-permeability boundary condition (U (t, x) · nD(x)) = 0

       Xt = X0 + t

0 Us ds,

Ut = U0 + t

0 B[Xs; ρ(s)]ds + σWt−

0<s≤t

2 (Us− · nD(Xs)) nD(Xs)1{Xs∈∂D}, ρ(t) is the probability density of (Xt, Ut) for all t ∈ (0, T], QT = (0, T) × D × Rd and ΣT = (0, T) × ∂D × Rd

Definition - Trace of the density along ΣT

γ(ρ) : ΣT → R is the trace of (ρ(t); t ∈ [0, T]) along ΣT if it is nonnegative and satisfies, for all t in (0, T], f in C∞

c

Qt
:
Σt

(u · nD(x)) γ(ρ)(s, x, u)f(s, x, u) ds dσ∂D(x) du =

D×Rd (f(0)ρ(0) − f(t)ρ(t)) +
Qt
∂sf + u · ∇xf + B[·; ρ·] · ∇uf + σ2

2 △uf

ρ(s)

and, for dt ⊗ dσ∂D a.e. (t, x) in (0, T) × ∂D,

Rd |(u · nD(x))|γ(ρ)(t, x, u) du < +∞,
Rd γ(ρ)(t, x, u) du > 0.

71

SLIDE 72

Theorem [B. and Jabir 14, 15]

Under (H), there exists a unique solution in law to the SLM with jumps in M(C([0, T]; D) × D([0, T]; Rd)), and this law admits time marginal densities ρt = P ◦ (x(t), u(t))−1 in L2(ω; D × Rd), for all t ∈ [0, T], that solves in V1(ω; QT)          ∂tρ(t, x, u) + u · ∇xρ(t, x, u) − σ2 2 △uρ(t, x, u) = − (B[x; ρ(t)] · ∇uρ(t, x, u)) in QT ρ(0, x, u) = ρ0(x, u) in D × Rd, γ−(ρ)(t, x, u) = γ+(ρ)(t, x, u − 2(u · nD(x))nD(x)) in Σ−

T ,

The trace γ(ρ) (in the sense of the previous definition) satisfies the no-permeability boundary condition E{(Ut · nD(Xt))|Xt = x} =

Rd(u · nD(x))γ(ρ)(t, x, u) du
Rd γ(ρ)(t, x, u) du

= 0, dt ⊗ dσ∂D

a.e. on (0, T) × ∂D.

The associated smoothed N-particles system propagates chaos.

72

SLIDE 73

The hypotheses (H)

(HLangevin) for the construction of the linear Langevin process (HMVFP) for the well-posedness of the Vlasov-Fokker-Planck equation

◮ (HLangevin)-(i) (X0, U0) is distributed according to the initial law µ0 having its

support in D × Rd and such that

D×Rd
|x|2 + |u|2

µ0(dx, du) < +∞.

◮ (HLangevin)-(ii) The boundary ∂D is a compact C3 submanifold of Rd. ◮ (HMVFP)-(i) b : Rd → Rd is a bounded continuous function. ◮ (HMVFP)-(ii) The initial law µ0 has a density ρ0 in the weighted space

L2(ω, D × Rd) with ω(u) := (1 + |u|2)

α 2

for some α > d + 3.

◮ (HMVFP)-(iii) There exist two measurable functions P0, P0 : R+ −

→ R+ such that 0 < P0(|u|) ≤ ρ0(x, u) ≤ P0(|u|), a.e. on D × Rd; and

Rd(1 + |u|)ω(u)P

2 0(|u|)du < +∞.

73

SLIDE 74

Propagation of chaos

With independent copies {(Xi

0, Ui 0, (Wi)), i = 1, . . . , N} of (X0, U0, (W))

                 Xi,ǫ,N

t

= Xi

0 +

t Ui,ǫ,N

s

ds, Ui,ǫ,N

t

= Ui

0 +

t Bǫ[Xi

s; ¯

µǫ,N

s

]ds + σWi

t + Ki,ǫ,N t

, Ki,ǫ,N

t

= −2

0<s≤t
Ui,ǫ,N

s−

· nD(Xi,ǫ,N

s

)

nD(Xi,ǫ,N

s

)1 Xi,ǫ,N

s

∈ ∂D , i = 1, . . . , N. where ¯ µǫ,N

t

= 1

N N

i=1

δ(Xi,ǫ,N

t

,Ui,ǫ,N

t

), Bǫ[x; γ] =

D×Rd b(v)βǫ(y)φǫ(x − y)γ(dy, dv)
D×Rd βǫ(y)φǫ(x − y)γ(dy, dv) + ǫ

Theorem

Let P be the law on E of (X, U, K), and let Pǫ,N be the law of {(Xi,ǫ,N, Ui,ǫ,N, Ki,ǫ,N); 1 ≤ i ≤ N}. Then Pǫ,N is P-chaotic; namely, for all {Fl, 1 ≤ l ≤ k}, k ≥ 2, with Fi ∈ Cb(C([0, T]; D) × D([0, T]; Rd) × D([0, T]; Rd)), lim

ǫ→0+

lim

N→+∞F1 ⊗ F2 ⊗ · · · ⊗ Fk ⊗ 1 ⊗ 1 ⊗ · · · , Pǫ,N = k

l=1

Fl, P.

74

SLIDE 75

Incompressible stochastic Lagrangian model in the torus

B., Fontbona, Jabin and Jabir 2013

Goal : well-posedness for the following process in Td × Rd: Xt = [X0 + t Usds] modulo 1 Ut = U0 + t

β (Us − αUs) − ∇P(s, Xs)

ρ(s, Xs) 1{ρ(s,Xs)>0}

ds + σWt.

with Law(Xt) = ρ(t, x)dx, Ut = E(Ut|Xt) and − △xP(t, x) =

d

i,j=1

∂2

ij

E(Ui

tUj t|Xt = x)ρ(t, x)

,

(t, x) ∈ R+ × Td.

75

SLIDE 76

Incompressible stochastic Lagrangian model in the torus

Lemma

Assume that the previous nonlinear SDE has a solution (X, U) such that E|Ut|2 < ∞ ∀t ∈ [0, T], E T |Us|2ds < ∞ and T

Td |∇P(s, x)|dxds < ∞.

Assume also that for all 1−periodic function ϕ : Rd → R of class C1 we have E(∇ϕ(X0) · U0) = 0. Then, a) E(∇ϕ(Xt) · Ut) = 0 for all t ∈ [0, T], b) the process (Xt, t ∈ [0, T]) is stationary.

Existence result for the Vlasov Fokker Planck equation ?

76

SLIDE 77

Local existence of analytical solutions

                     ∂tf(t, x, u) + u · ∇xf(t, x, u) = σ2 2 △uf(t, x, u) + β d f(t, x, u) + βu · ∇uf(t, x, u) + ∇uf(t, x, u) ·

∇P(t, x) − βα
Rd vf(t, x, v) dv
= 0 on (0, T] × Td × Rd,

f(t = 0) = f0 on Td × Rd, △P(t, x) = −∂xixj

Rd vivjf(t, x, v)dv, on [0, T] × Td,

Theorem (d=1)

Let ¯ λ > 0 and s ≥ d + 2 be an even integer. There exists κ0 = κ0(¯ λ, s) and r → κ1(r, ¯ λ, s) ≥ 0 s.t. if f0 : Td × Rd → R of class C∞ and T > 0 satisfy:

◮ Rd f0(x, u)du = 1 and ∇x ·

Rd uf0(x, u)du = 0 for all x ∈ T,

◮ (1 + |u|2)

s 2 Dl

xDk uf0∞ ≤ C0(|k|+m)!(|l|+n)! λ

|k|+|l|

for some n, m ∈ N, all pair of multiindexes k, l ∈ NN and a constant C0 < κ0(¯ λ, s), and

◮ T < κ1(C0, ¯

λ, s), then, a solution f of class C1,∞ exists.

77

SLIDE 78

Outline Stochastic particle methods for McKean SDEs Lagrangian viewpoint in turbulence modeling

Applications Atmospheric boundary layer simulation Wind farm simulation

Well posedness results of SLM (toy models)

boundary condition mass constraint

Numerical analysis for SLM (toy models)

78

SLIDE 79

We simplify again the fluid particle model Goal : analyze the rate of convergence

       Xt = X0 + t Us ds Ut = U0 + t B[Xs; ρs] ds + Wt, ρt = L(Xt, Ut), for all t ≤ T Rd × L1(Rd × Rd) ∋ (x, γ) → B[x; γ] =

Rd b(v)γ(x, v) dv
Rd γ(x, v)dv

1{

Rd γ(x,v)dv>0}

(t, x) → B[x; ρ(t)] = E[b(Ut)

Xt = x]

The smoothed process (Xǫ, Uǫ, ǫ > 0) converges in law to (X, U), by substituting the conditional expectation B[x; ρt] with the kernel regression estimate [Bossy Jabir Talay 2011] Bǫ[x; ρ] :=

Rd×Rd b(v)Kǫ(x − y) ρ(dy, dv)
Rd×Rd Kǫ(x − y) ρ(dy, dv)

, where Kǫ(x) := 1 ǫd K(x ǫ) and K is a kernel function.        Xǫ

t = X0 +

t Uǫ

s ds

Uǫ

t = U0 +

t Bǫ[Xǫ

s , ρǫ s] ds + Wt,

ρǫ

t = L(Xǫ t , Uǫ t ),

for all t ≤ T.

79

SLIDE 80

We simplify the fluid particle model Goal : analyze the rate of convergence

       Xt =

X0 +

t Us ds

mod 1

Ut = U0 + t B[Xs; ρs] ds + Wt, ρt = L(Xt, Ut), for all t ≤ T T × L1(T × Rd) ∋ (x, γ) → B[x; γ] =

Rd b(v)γ(x, v) dv
Rd γ(x, v)dv

1{

Rd γ(x,v)dv>0}

(t, x) → B[x; ρ(t)] = E[b(Ut)

Xt = x]

The smoothed process (Xǫ, Uǫ, ǫ > 0) converges in law to (X, U), by substituting the conditional expectation B[x; ρt] with the kernel regression estimate Bǫ[x; ρ] :=

T×Rd b(v)Kǫ(x − y) ρ(dy, dv)
T×Rd Kǫ(x − y) ρ(dy, dv)

, where Kǫ(x) := 1 ǫd K(x ǫ) and K is a kernel function.        Xǫ

t =

X0 +

t Uǫ

s ds

mod 1

Uǫ

t = U0 +

t Bǫ[Xǫ

s , ρǫ s] ds + Wt,

ρǫ

t = L(Xǫ t , Uǫ t ),

for all t ≤ T.

80

SLIDE 81

Numerical analysis of a non degenerate toy model (d=1)

Yt = Y0 + t b[Ys; ρs]ds + σwt, ρt = L(Yt). Assume b symmetric b(x, y) = b(x − y). b[x, ρt] = Eb(x − Yt). Introduce N particles: Let (wi, i = 1, . . . , N) N-Brownian motion, independent of the i.i.d. (Yi

0, i = 1, . . . , N) with

law ρ0. We consider Yi,N

t

= Yi

0 +

t 1 N

N

j=1

b(Yi,N

s

− Yj,N

s )ds + σwi t,

µN

t = 1

N

i=1

δYi,N

t

.

Lemma

Assume that there exists r ≥ 2 such that b(·) is in Cr+1

b

(R), and ρ0 ∈ Wr,∞(R) ∩ Wr,1(R). Then ρt ∈ Wr,∞, and for all test function φ ∈ C3

c(R), there exists a positive constant C

such that sup

t∈[0,T]

E

φ, µN

t − ρt

≤

C √ N . where the constant C is of the form β exp(αb′ t).

81

SLIDE 82

Control of the drift error term

∀t ∈ [0, T], Edrift(t) := E

b[Y1,N

t

, ρt] − 1 N

N

j=1

b(Y1,N

t

− Yj,N

t

)

.

Notice that, with Yi

t = Yi 0 +

t

0 b[Yi s; ρs]ds + σwi t,

E|Y1

t − Y1,N t

| =E

t

 b[Y1

s , ρs] − 1

N

j=1

b(Y1,N

s

− Yj,N

s

)   ds

Since b is Lipschitz with constant L = b′∞, x → b[x, ρt] is Lipschitz, and applying Gronwall

Lemma E|Y1

t − Y1,N t

| ≤ t LE|Y1

s − Y1,N s

|ds + t Edrift(s)ds ≤ t exp(L(t − s))Edrift(s)ds. Edrift(t) ≤ 3LE|Y1,N

t

− Y1

t | + E

b[Y1

t , ρt] − 1

N

j=1

b(Y1

t − Yj t)

≤

t 3L exp(L(t − s))Edrift(s)ds + E

Eb(x − Yt) − 1

N

j=1

b(x − Yj

t)

x=Y1

t

.

The (Yj, j = 1) being i.i.d. E

Eb(x − Yt) − 1

N N

j=1 b(x − Yj t)

2

x=Y1

t

≤ C

N . Thus Edrift(t) ≤ t 3L exp(L(t − s))Edrift(s)ds + C √ N ≤ exp(4Lt) C exp(−Lt) √ N ≤ exp(3Lt) C √ N .

82

SLIDE 83

A symptomatic case

Yt = Y0 + t 1 2ρs(Ys)ds + wt, ρt = L(Yt).

McKean 1967, Calderoni &Pulvirenti 1983, Sznitman & Varadhan 1986, Sznitman 1986.

Yi,N

t

= Yi

0 +

t Kǫ[Yi,N

s ; µN s ]ds + σwi t,

µN

t = N

i=1

δYi,N

t

with Kǫ[x; γ] :=

R Kǫ(x − y)γ(dy) with Kǫ(z) = 1

ǫK( z ǫ) and K positive function of mass

equal to 1. When ρ0 ∈ W2,1(R) ∩ W2,∞(R) then ρt ∈ W2,1(R) ∩ W2,∞(R) for all finite t > 0 and ρt − ρǫ

t L1(R) ≤ Cǫ2,

where Yǫ

t = Y0 +

t

1 2Kǫ[Yǫ s ; ρǫ s]ds + wt,

ρǫ

t = L(Yǫ t ).

From the previous lemma: supt∈[0,T] E

φ, µN,ǫ − ρt
≤ C
ǫ2 + exp(C/ǫ)

√ N

.

Numerical experiments are more optimistic :

0.1 0.2 0.3 0.4 0.5 0.6 0.7

2.5
2
1.5
1
0.5

0.5 1 1.5 2 2.5 3 BURGERS : pas=0.1 / 16000 particules/ T= 0.2 /Diffusion= 1 / Epsilon = 0.500E-2 solution exacte Approximation 0.1 0.2 0.3 0.4 0.5 0.6

2.5
2
1.5
1
0.5

0.5 1 1.5 2 2.5 3 BURGERS : pas=0.1 / 16000 particules/ T= 0.2 /Diffusion= 1 / Epsilon = 0.100 solution exacte Approximation

(ρ0(dx) = 1[0,1] with N = 16000, T = 0.2, σ = 1, ǫ = 0.510−2 left, ǫ = 0.1 right)

83

SLIDE 84

Go back to McKean kinetic particle system

(Xi

t, Ui t, t ≤ T, i = 1, . . . , N) are the strong solution of

       Xi

t =

Xi

0 +

t Ui

s ds

mod 1

Ui

t = Ui 0 +

t B[Xi

s, ρs] ds + Wi t,

ρt = L(Xi

t, Ui t),

for all t ≤ T, where we recall that (Xi

0, Ui 0, i = 1, . . . , N) are i.i.d. with density ρ0.

(Xi,ǫ

t , Ui,ǫ t , t ≤ T, i = 1, . . . , N) are the strong solution of

       Xi,ǫ

t

=

Xi

0 +

t Ui,ǫ

s ds

mod 1

Ui,ǫ

t

= Ui

0 +

t Bǫ[Xi,ǫ

s , ρǫ s] ds + Wi t,

ρǫ

t = L(Xi,ǫ t , Ui,ǫ t ),

for all t ≤ T,                  Xi,N

t

= Xi

0 +

t Ui,N

s ds,

Ui,N

t

= Ui

0 +

t

1 N

N

j=1 b(Uj,N s )Kǫ(Xi,N s

− Xj,N

s ) 1 N

N

j=1 Kǫ(Xi,N − Xj,N s )

Nadaraya-Watson estimator at Xi,N

s

ds + Wi

t,

for all t ≤ T. Control the error seen by the particles : E

F[X1

t ; ρt] − Fǫ[X1,N t

; µN

t ]

.

84

SLIDE 85

Rate of convergence

For any f measurable bounded on Rd, we set the kernel regression version of the conditional expectation : (x, t) →F[x; ρt] = E[f(Ut)|Xt = x], (x, t) →Fǫ[x; ρt] =

Rd×Rd f(v)Kǫ(x − y) ρt(dy, dv)
Rd×Rd Kǫ(x − y) ρt(dy, dv)

= E[f(Ut)Kǫ(x − Xt)] E[Kǫ(x − Xt)] Assume that f and b smooth and bounded function with bounded derivatives; K > 0 Lipschitz bounded, with compact support.

Theorem

For any 1 < p < 1 + 1 1 + 3d and c > 0, there exists a constant C such that for all α > 0, ǫ > 0 and N > 1 satisfying 1 α2ǫdN

1 p

≤ c, we bound the error seen by the particles by E

F[X1

t ; ρt] − Fǫ[X1,N t

; µN

t ]

≤ C
ǫ +

1 αǫdN + 1 ǫ(d+1)pN + 1 (ǫdN)

1 p

+ 1 ǫ

dp 2 √

N

.

If we choose α = ǫ, the optimal rate of convergence is achieved for N = ǫ−(d+2)p and E

F[X1

t ; ρt] − Fǫ[X1,N t

; µN

t ]

≤ N

−

1 (d+2)p . 85

SLIDE 86

Error decomposition Smoothing error : there exists C, such that for all t ∈ [0, T] for all ǫ > 0.

E|F[Xt; ρt] − Fǫ[Xǫ

t ; ρǫ t ]| ≤ Cǫ

Kernel regression estimates error : there exists a constant C such that for

all t ≤ T, ǫ > 0 and N > 0,

T

E  

Fǫ[z; ρǫ

t ] −

j=i f(Uj,ǫ

t )Kǫ(z − Xj,ǫ t )

j=i Kǫ(z − Xj,ǫ

t ) + (N − 1)α

2

 ρǫ

t (z) dz ≤ C

1

ǫd(N − 1) + ǫ2

.

Statistical error : for any p > 1 sufficiently close to 1, there exists a

constant C such that for α > 0, ǫ > 0, N > 1 and t ≤ T E

Fǫ[Xi,N

t

; ρǫ

t ] −

N

j=1 Eǫ 0,Xj

0,Uj 0[f(Uǫ

t )Kǫ(Xi,N t

− Xǫ

t )]

N

j=1 Eǫ 0,Xj

0,Uj 0[Kǫ(Xi,N

t

− Xǫ

t )] + Nα

≤ Cǫ +

C ǫd+1N + C αǫdN .

The iteration error

86

SLIDE 87

The iteration error

Step 1 : remove the denominators E

N

j=1 Eǫ 0,Xj

0,Uj 0[f(Uǫ

t )Kǫ(Xi,N t

− Xǫ

t )]

N

j=1 Eǫ 0,Xj

0,Uj 0[Kǫ(Xi,N

t

− Xǫ

t )] + Nα

− N

j=1 f(Uj,N t

)Kǫ(Xi,N

t

− Xj,N

t

) N

j=1 Kǫ(Xi,N t

− Xj,N

t

) + Nα

≤C

N E

N
j=1
Eǫ

0,Xj

0,Uj 0[f(Uǫ

t )Kǫ(Xi,N t

− Xǫ

t )] − f(Uj,N t

)Kǫ(Xi,N

t

− Xj,N

t

)

+ C

N E

N
j=1
Eǫ

0,Xj

0,Uj 0[Kǫ(Xi,N

t

− Xǫ

t )] − Kǫ(Xi,N t

− Xj,N

t

)

+

C (ǫdN)

1 p

Step 2 : control the derivative of u(k)

t,χ(s, y, v) = Es,y,v[b(k)(Uǫ t )Kǫ(χ − Xǫ t )].

For all t ≤ T, and p such that 1 ≤ p < 1 +

1 1+3d, we have

sup

x,u

dχ |∇vuT,χ(t, x, u)|p

p ≤

CT (T − t)

p 2 + 3 2 d(p−1) .

for all ǫ > 0. In particular, this quantity is time integrable on [0, T] and the bound of the time integral is independent of ǫ.

87

SLIDE 88

The iteration error

Step 3 : particle removal in 1 N E

j=i

t ∇vu(k)

t,Xi,N

t (s, Xj,N

s , Uj,N s ) ·

Bǫ[Xj,N

s ; µN s ] − Bǫ[Xj,N s ; ρǫ s]

ds
,

Let p > 1 be sufficiently small. There exists a constant C such that for all α, ǫ and N satisfying ǫdN > 1, t E

|∇vu(k)

t,Xi,N

t (s, Xj,N

s , Uj,N s )|1 ·

B(l)

ǫ [Xj,N s ; µN s ] − B(l) ǫ [Xj,N s ; ρǫ s]

ds

≤ t C (t − s)

p 2 E

|B(l)

ǫ [Xj,N−1 s

; ρǫ

s] − B(l) ǫ [Xj,N−1 s

; µN−1

s

]|

ds

+ C α2(ǫdN)

1 p

t E

1

N − 1

N−1
j=1

Kǫ(Xk,N−1

s

− Xj,N−1

s

) − Eǫ

0,Xj

0,Uj 0[Kǫ(Xk,N−1

s

− Xǫ

s )]

ds

+ C (ǫdN)

1 p

1 +

1 α2(ǫdN)

1 p

+

C αǫdN , for all t ≤ T and all j = i.

88

SLIDE 89

The iteration error

Step 4 : particle inclusion For any p > 1 and c > 0, there exists a constant C such that for all s and t satisfying 0 ≤ s < t ≤ T and all α > 0, ǫ > 0, N > 1 satisfying 1 α2(ǫdN)

1 p

≤ c, we have E

1

N−1

N−1
j=1

Kǫ(Xk,N−1

t

− Xj,N−1

t

) − Eǫ

0,Xj

0,Uj 0[Kǫ(Xk,N−1

t

− Xǫ

t )]

≤ C

ǫ

d q

sup

s≤t

E

1

N−1

N−1
j=1

Kǫ(Xk,N

s

− Xj,N

s ) − Eǫ 0,Xj

0,Uj 0[Kǫ(Xk,N

s

− Xǫ

s )]

+

C (ǫdN)

1 p

E

|Bǫ[Xj,N−1

t

; ρǫ

t ] − Bǫ[Xj,N−1 t

; µN−1

t

]|

≤ CE
|Bǫ[Xj,N

t

; ρǫ

t ] − Bǫ[Xj,N t

; µN−1,N

t

]|

+ C sup

s≤t

E

1

N−1

N−1
j=1

Kǫ(Xk,N

s

− Xj,N

s ) − Eǫ 0,Xj

0,Uj 0[Kǫ(Xk,N

s

− Xǫ

s )]

+

C (ǫdN)

1 p

+ C αǫdN , where µN−1,N is the normalized empirical distribution of the N − 1 processes (Xk,N, Uk,N) for all k ≤ N − 1: µN−1,N =

1 N−1

N−1

j=1 δ{(Xj,N,Uj,N)}.

89

SLIDE 90

Smoothing error

E

F[Xt; ρt] − Fǫ[Xǫ

t ; ρǫ t ]

≤ E|F[Xt; ρt] − Fǫ[Xt; ρǫ

t ]|

A

+ E|Fǫ[Xt; ρǫ

t ] − Fǫ[Xǫ t ; ρǫ t ]|

B

. A: use the uniform density lower bound + kinetic semigroup properties B: use moreover that the kernel function K is in C1

b(Rd) and there exists a

constant C such that |∇Kǫ(x)| ≤ C ǫ Kǫ(x) for all x and ǫ > 0. Moreover, assume there exists a constant αT such that for all ǫ > 0 and t ≤ T, we have αT ≤ inf

χ ρǫ t ⋆ Kǫ(χ),

Then, there exists a constant C such that for all ǫ > 0, sup

t≤T

sup

x∈T

|∇Fǫ[x; ρǫ

t ]| ≤ C.

90

SLIDE 91

Rate of convergence for the measured error

For any f measurable bounded on Rd, we set the kernel regression version of the drift Fǫ[x; ρ] :=

Rd×Rd f(v)Kǫ(x − y) ρ(dy, dv)
Rd×Rd Kǫ(x − y) ρ(dy, dv)

, (x, t) → F[x; ρt] = E[f(Ut)|Xt = x] f and b smooth and bounded function with bounded derivatives

Proposition

For any 1 < p < 1 +

1 1+3d and c > 0, there exists a constant C such that for all α > 0,

ǫ > 0 and N > 1 satisfying 1 α2ǫdN

1 p

≤ c, we have E|Fǫ[Xj,N

t

; ρǫ

t ] − Fǫ[Xj,N t

; µN

t ]

≤ C
ǫ +

1 αǫdN + 1 ǫ(d+1)pN + 1 (ǫdN)

1 p

+ 1 ǫ

dp 2 √

N

.

If we choose α = ǫ, the optimal rate of convergence is achieved for N = ǫ−(d+2)p and E

Fǫ[Xj,N

t

; ρǫ

t ] − Fǫ[Xj,N t

; µN

t ]

≤ N

−

1 (d+2)p . 91

SLIDE 92

Regression Kernels

We consider only boxed kernel K such that b1S0,r(x) ≤ K(x) ≤ 1S0,R

92

SLIDE 93

Nonparametric regression estimate

Let (Xi, Ui) be an i.i.d sequence of n random variables.

Local averaging estimate :

approximate F(x) = E[f(U)|X = x] with: Fn(x) =

n

i=1

Wn,i(x) · f(Ui), where the weight Wn,i depends only on x and (Xj, f(Uj)), j ≤ n.

Kernel estimate :

Wn,i(x) = Kǫ(x − Xi) n

j=1 Kǫ(x − Xj),

where Kǫ(x) = K( x

ǫ). Related to the minimization:

Fn(x) = arg min

c

1 n

n

i=1

Kǫ(x − Xi)(f(Ui) − c)2 .

93

SLIDE 94

Nonparametric regression estimate

Partitioning estimate :

Let Pn = {An,1, An,2, . . . } be a partition of the domain. Set Wn,i(x) = 1{Xi∈An,j} n

i=1 1{Xi∈An,j}

, for x ∈ An,j. Related to the least squares estimate: Fn = arg min

fn

1 n

n

i=1

(f(Ui) − fn(Xi))2 , where the arg min is taken on the set of piecewise constant functions adapted to Pn. Also: local modeling estimates (e.g. local polynomial estimates), global modeling/least squares estimates, penalized modeling. . .

94

SLIDE 95

Convergence rate of the nonparametric estimates

◮ Consider the mean L2 error:

E

D

F(x) − Fn(x)
2ρ(x) dx
.

◮ Variance-Bias decomposition:

E

D

F(x) − Fn(x)
2ρ(x) dx
=
Var Fn(x)ρ(x) dx +
|bias(x)|2ρ(x) dx.

◮ The previous estimators are weakly universally consistent (Stone’s

theorem).

◮ For (p, C) smooth conditional expectations, the lower optimal convergence

rate is of order n−

2p 2p+d .

Theorem

For Lipschitz continuous regression functions, we have EFn − F2 ≤ C 1 nǫd + dCǫ2, for both partitioning estimate and kernel estimate (with naive kernel). This gives an optimal rate of n−

2 d+2 . 95

SLIDE 96

Local averaging based particle algorithm

Particle algorithm :

1. Initialization of (Xi

0, Ui 0) for 1 ≤ i ≤ n

2. Time loop : for tk = k∆t up to T :

◮ for each particle 1 ≤ i ≤ n :

2.1 update the position : Xi,n

k

= Xi,n

k−1 + Ui,n k−1

2.2 add gaussian noise: Ui,n

k

= Ui,n

k−1 + (Wi k − Wi k−1)

2.3 drift computation : loop over each particle to calculate the Wn,j(Xi,n

k−1):

Ui,n

k + =

N

j=1 b(Uj,n k−1)Kn(Xi,n k−1 − Xj,n k−1)

N

j=1 Kn(Xi,n k−1 − Xj,n k−1)

∆t

◮ End loop particles

3. Final mean field evaluation at x: Loop on all particles

n

j=1 f(Uj,n T )Kn(x − Xj,n T )

n

j=1 Kn(x − Xj,n T )

Complexity up to O(n2) for kernel estimates, and only O(n) for partitioning estimates.

96

SLIDE 97

Numerical simulation

For d = 2

Test case

dUt = Cdt − ∇V(Xt)dt

“pressure”

+ (Ut − 2Ut)dt

conditional drift

+ dWt, with V(x, y) = 1 2π cos(2π x) sin(2π y) − 1 2x, X0 = σZ with L(Z) = N(0, 1), L(U0) = N(0, 1), σ = 0.3.

Reference simulation

Regression estimate implemented with Epanechnikov kernel T = 2, 128 times step, ¯ N = 105 particles, ǫ =

1 16.

97

SLIDE 98

98

SLIDE 99

We want to observe the error of the scheme as

D

|F[x; ρT] − Fǫ[x; µǫ,N

T ]|ρT(x) dx,

First we approximate the mean field F[x; ρT] by Fǫ

x; µǫ,¯

N] :=

1 Nmc

Nmc

k=1

Fǫ

x; µǫ,¯

N(ωk)

,

Using spline interpolation on a mesh of size ∆x, we observe the approximate L1error 1 Nmc

Nmc

k=1
D
Fǫ
x; µǫ,¯

N] ∆x

− F∆x

ǫ [x; µǫ,N(ωk)]

¯

ρ∆x(x) dx.

99

SLIDE 100

Observe the approximate L1error:

10−1 ǫ 105 N . 1 3 0.014 0.015 0.015 0.016 0.017 0.018 0.019 0.020 0.021 0.022 0.022 0.023 0.023 0.024 0.024 . 2 5 0.025 0.026 0.026 0.028 0.029 . 3 1

Optimal N(ǫ) N = 1/ǫ4

in-line with what we expected from the theoretical rate where the optimal value

f window size is given by

1 N

1 d+2 . 100

SLIDE 101

Observe the kernel regression estimation error

bias :

D
Fǫ
x; µǫ,N∆x

− Fǫ

x; µǫ,¯

N∆x

2

¯ ρ∆x

ǫ (x) dx,

10−1 ǫ 10−2 MDIB

N = 10000 N = 14962 N = 22388 N = 33499 N = 50124 N = 75000 ǫ−1.0 ǫ−1.0

variance : 1 Nmc

Nmc

j=1
D
Fǫ
x; µǫ,N∆x

− F∆x

ǫ

x; µǫ,N(ωj)
2

¯ ρ∆x

ǫ (x) dx.

10−1 ǫ 10−2 10−1 MDIVar

N = 10000 N = 14962 N = 22388 N = 33499 N = 50124 N = 75000 ǫ−1.0 ǫ−1.0

101

SLIDE 102

Impact of the kernels

Velocity field and its error

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.40 0.48 0.56 0.64 0.72 0.80 0.88 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 −0.016 −0.008 0.000 0.008 0.016 0.024 0.032 0.040 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.40 0.48 0.56 0.64 0.72 0.80 0.88 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 −0.016 −0.008 0.000 0.008 0.016 0.024 0.032 0.040 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.40 0.48 0.56 0.64 0.72 0.80 0.88 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 −0.016 −0.008 0.000 0.008 0.016 0.024 0.032

Regression kernels : hat kernel triangle kernel epanechnikov kernel

102

SLIDE 103

Impact of the kernels

Partitioning kernels :

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.36 0.42 0.48 0.54 0.60 0.66 0.72 0.78 0.84 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 −0.16 −0.12 −0.08 −0.04 0.00 0.04 0.08 0.12 0.16

Indicator function (NGP)

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 0.40 0.48 0.56 0.64 0.72 0.80 0.88 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 −0.02 −0.01 0.00 0.01 0.02 0.03 0.04 0.05

Hat function (CIC)

103

SLIDE 104

Coauthors

On stochastic analysis of SLM, numerical analysis Jean Francois Jabir (University of Valpareiso ) Joaquin Fontbona (U Chile) Pierre Emmanuel Jabin (U Maryland) Laurent Violeau (Former Phd student Inria) SDM and WindPoS projects http://windpos.inria.fr Antoine Rousseau (Inria Lemon) Philippe Drobinski (LMD-IPSL) Jos´ e Espina, Jacques Morice, Cristian Paris (Inria Chile) Selim Kraria (Inria SED)

Thank you

SLIDE 105

Bibliography

F . Bernardin, M. Bossy, C. Chauvin, P . Drobinski, A. Rousseau, and T. Salameh. Stochastic Downscaling Methods : Application to Wind Refinement. Stoch. Environ. Res. Risk. Assess., 23(6), 2009. F . Bernardin, M. Bossy, C. Chauvin, J-F. Jabir, and A. Rousseau. Stochastic Lagrangian Method for Downscaling Problems in Computational Fluid Dynamics. ESAIM: M2AN, 44(5):885–920, 2010.

M. Bossy, J. Fontbona, P-E. Jabin, and J-F. Jabir.

Local existence of analytical solutions to an incompressible lagrangian stochastic model in a periodic

domain. Communications in Partial Differential Equations, 38(7):1141–1182, 2013.
M. Bossy, J.F

. Jabir, and D. Talay. On conditional McKean Lagrangian stochastic models. Probab. Theory Relat. Fields, 151:319–351, 2011.

M. Bossy and J-F. Jabir.

On confined mckean langevin processes satisfying the mean no-permeability boundary condition. Stochastic Processes and their Applications, 121(12):2751 – 2775, 2011.

M. Bossy and J.F

. Jabir. Lagrangian stochastic models with specular boundary condition. Journal of Functional Analysis, 268(6):1309 – 1381, 2015.

L. Violeau. Particles method for Lagrangian Stochastic Models.

PhD thesis, PhD Thesis, Universit´ e de Nice Sophia Antipolis. In progress.

M. Bossy, J. Espina, J. Morice, C. Paris, and A. Rousseau.

Modeling the wind circulation around mills with a lagrangian stochastic approach. Preprint arXiv:1404.4282, 2014.