Stochastic Lagrangian Models Mireille Bossy Tosca Laboratory - PowerPoint PPT Presentation

The tow following assertions are equivalent : (Tanaka 82) • P N is P -chaotic • µ N converges weakly toward the (deterministic) measure P Theorem (M´ el´ eard 95) � � � � | X i , N | 2 | X i t | 2 sup E sup + sup E sup < ∞ t N 0 ≤ t ≤ T N 0 ≤ t ≤ T � � | X i , N − X i t | 2 < ∞ , sup N E sup t N 0 ≤ t ≤ T where � t � t X i t = X i σ [ X i s , ρ s ] dW i b [ X i 0 + s + s , ρ s ] ds 0 0 with ρ t = P ◦ X − 1 . t 17

Numerical approximation � t � t X t = X 0 + σ [ X s , ρ s ] dW s + b [ X s , ρ s ] ds , with ρ s = Law ( X s ) . 0 0 Spatial discretisation : � t � t X i , N X i σ [ X i , N s , µ N s ] dW i b [ X i , N s , µ N s ] ds , i ≤ N . = 0 + s + t 0 ) Time-Discretisation : ∆ t > 0 , such that T = K ∆ t for K ∈ N ,  � N   k ∆ t + ( 1 i , N j , N  X i , N X i , N ¯ ( k + 1 )∆ t = ¯ k ∆ t ))( W i ( k + 1 )∆ t − W i  σ ( X k ∆ t , X k ∆ t )  N    j = 1 N � + ∆ t ( 1  b (¯ X i , N k ∆ t , ¯ X j , N k ∆ t )) ,    N   j = 1   X i , N ¯ t = 0 = X i 0 . � N ( x ) = 1 N ,ε, ∆ t φ ε ( x − ¯ X i , N k ∆ t ) , pour un cut-off donn´ e φ . U k ∆ t N i = 1 18

Rate of convergence Theorem (Bossy & Talay 97) Dimension d = 1 . K ( R ) , b ( · , · ) ∈ C 2 , 2 Assume that ρ 0 ( · ) ∈ C 2 b ( R 2 ) and s ( · , · ) ∈ C 3 b ( R 2 ) . Assume the strong ellipticity of the differential operator L . Then ρ t ( · ) ∈ W 2 , 1 ( R ) for all t ∈ [ 0 , T ] . For a 2-order cut-off φ ( x ) , ∀ k ≤ K , � � √ ∆ t 1 N ,ε, ∆ t N + ε 2 + √ E � U − ρ k ∆ t � L 1 ( R ) ≤ C b ,σ, T ,ρ 0 . k ∆ t ε ε H ( x ) := 1 { x ≥ 0 } , Heaviside function � 1 � N � √ E � ( H ∗ ρ k ∆ t )( x ) − 1 i , N H ( x − X k ∆ t ) � L 1 ( R ) ≤ C b ,σ, T ,ρ 0 √ N + ∆ t . N i = 1 19

Rate of convergence Theorem (Antonelli & Kohatsu-Higa 02) d = 1 . b ( · , · ) and s ( · , · ) C ∞ ( R 2 ) , with bounded derivatives. A non degenerescency condition for L (Malliavin calculus). φ ( x ) is the Gaussian density on R . ε 2 = ∆ t . � � 1 1 N ,ε, ∆ t √ √ E � U − ρ k ∆ t � L 1 ( R ) ≤ C b ,σ, T ,ρ 0 + ∆ t + . k ∆ t ∆ tN N H ( x ) := 1 x ≥ 0 , Heaviside function, � 1 � � N E � ( H ∗ ρ k ∆ t )( x ) − 1 i , N H ( x − X k ∆ t ) � L 1 ( R ) ≤ C b ,σ, T ,ρ 0 √ N + ∆ t . N i = 1 20

Incompressible 2 d -Navier Stokes equation V ( t , x ) = ( v 1 ( t , x ) , v 2 ( t , x )) and p ( t , x ) solution to  ∂ t V = ν ∆ V − ( V · ∇ ) V − ∇ p , t ≥ 0 , x ∈ R 2 ,  ∇ · V = ∂ x 1 v 1 + ∂ x 2 v 2 = 0 , (7)  V ( 0 , x ) = V 0 ( x ) . Vorticity formulation : U = ∇ × V = ∂ x 1 v 2 − ∂ x 2 v 1 obtained by derivation.  � � ∂ U  U ( � , x ∈ R 2 ,  ∂ t ( t , x ) = ν ∆ U ( t , x ) − ∇ . K ∗ U )( t , x ) (8) V ( x , t ) = ( � � K ∗ U )( x , t ) ,   U ( 0 , x ) = ( ∇ × V 0 ) ( x ) , 1 � with K ( x 1 , x 2 ) = 2 )( − x 2 , x 1 ) . 2 π ( x 2 1 + x 2 Chorin 73 : Random Vortex Method. 21

N. S. incompressible 2 d : probabilistic Interpretation Marchioro & Pulvirenti 82, Osada 85, M´ el´ eard 01 - Introduce a martingale problem, where the weighted marginal density laws ( � P t ) t ≥ 0 from a solution P , solve the vortex equation and V ( t , x ) = ( � K ∗ � P t )( 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 V N ( t , x ) = ( � µ N K ∗ � t )( x ) to V ( t , x ) . 22

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 U 0 has compact support. Fix h , Λ h = { h . i , i ∈ Z 2 } . � U h h 2 w i δ ( x − α i ) 0 ( x ) = Λ h ∩ Supp ( U 0 ) φ ∈ C L ( R 2 ) cut–off function of order m , with fast decay at infinity The particle system : � � � √ dX ε, h , i X ε, h , i − X ε, h , j w j h 2 K ε 2 ν dW i = dt + t . t t t j � � � V ε, h ( t , x ) = x − X ε, h , j w j h 2 K ε . t j 23

N. S. incompressible 2 d : rate of convergence ? (Long 88) Choose ε ≥ h q , q < 1 . Assume that V ( t , x ) is smooth enough. � � � h � L � � ε m + � V ε, h ( t , x )( ω ) − V ( t , x ) � sup L p ( B ( R 0 )) ≤ C ε + h | ln h | ε 0 ≤ t ≤ T except for ω in A , an event of probability smaller than h C ′ C , with C > 0 and C ′ > 0 depending in all the parameters ( T , L , m , p , q , R 0 , V 0 , ∂ i , j V ) . 24

Incompressible N.S. in a bounded domain O ⊂ R 2 V ( t , x ) = ( v 1 ( t , x ) , v 2 ( t , x )) and p ( t , x ) solve  ∂ t V = ν ∆ V − ( V · ∇ ) V − ∇ p , t ≥ 0 , x ∈ O ,  ∇ · V = ∂ x 1 v 1 + ∂ x 2 v 2 = 0 (9)  V ( 0 , x ) = V 0 ( x ) , V ( t , x ) = 0 , x ∈ ∂ O . Vorticity equation : U = ∂ x 1 v 2 − ∂ x 2 v 1  ∂ U  ∂ t ( t , x ) = ν ∆ U ( t , x ) − ∇ . ( UV ( t , x )) , x ∈ O ,  ∇ · V = 0 ⇒ V = ∇ ⊥ ψ = ( ∂ x 2 ψ, − ∂ x 1 ψ ) with U = − ∆ ψ,   U t = 0 = U 0 = ∇ × � V 0 , � 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 ) = − x γ ( x , y ) U ( t , y ) dy . O Resulting boundary condition for the vorticity : ∀ x ∈ ∂ O , ψ ( t , x ) = 0 ⇒ ∂ τ ψ = 0 ⇒ V . n = 0 . 25

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 , U t = 0 = U 0 = ∇ × � V 0 in O ,   ∂ n U ( 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

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

SLM for Atmospheric Boundary Layer simulation 28

Multiphase flows 29

Statistical approach of turbulent flows Eulerian description ◮ The incompressible Navier Stokes equation in R 3 , for the velocity field ( U ( 1 ) , U ( 2 ) , U ( 3 ) ) and the pressure P , with constant mass density ̺ ∂ t U ( ω ) + ( U ( ω ) · ∇ ) U ( ω ) = ν △ U ( ω ) − 1 ̺ ∇ P ( ω ) , t > 0 , x ∈ R 3 , ∇ · U ( ω ) = 0 , t ≥ 0 , x ∈ R 3 , U ( 0 , x ) ( ω ) = U 0 ( x ) ( ω ) , x ∈ R 3 . ◮ Ensemble averages as expectations: � � U � ( t , x ) := U ( t , x , ω ) d P ( ω ) . Ω 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

Reynolds-Averaged Navier-Stokes (RANS) Equation Assuming Reynolds decomposition, with constant mass density ̺ , we get � 3 � 3 j � = ν △� U ( i ) � − 1 ∂ t � U ( i ) � + � U ( j ) � ∂ x j � U ( i ) � + ∂ x j � u ′ i u ′ ̺∂ x i � P � , j = 1 j = 1 ∇ · � U � = 0 , t ≥ 0 , x ∈ R 3 , � U � ( 0 , x ) = � U 0 � ( x ) , x ∈ R 3 , closed by a chosen parameterization or a model equation of the Reynolds stress tensor ( � u ′ i u ′ j � = � U ( i ) U ( j ) � − � U ( i ) �� U ( j ) � , i , j ) . � 3 Turbulent kinetic energy tke ( t , x ) := 1 � u ′ i u ′ i � ( t , x ) 2 i = 1 � 3 � 3 � ∂ x k u ′ i ∂ x k u ′ i � ( t , x ) Pseudo-dissipation ε ( t , x ) := ν i = 1 k = 1 31

Equation for the Reynolds stress ( � u ′ i u ′ j � , i , j ) � 3 � � ∂ t � u ′ i u ′ � U � · ∇ x � u ′ i u ′ ∂ x k � u ′ i u ′ j u ′ j � + j � + k � k = 1 − 1 j ∂ x i p ′ + u ′ ̺ � u ′ i ∂ x j p ′ � + ν △ x � u ′ i u ′ = j � � �� � velocity pressure gradient tensor Π ij � 3 � ∂ x k u ′ i ∂ x k u ′ + ν j � k = 1 � �� � dissipation tensor ε ij � � 3 � � u ′ i u ′ k � ∂ x k � U ( i ) � + � u ′ j u ′ k � ∂ x k � U ( j ) � − k = 1 � �� � turbulence production tensor P ij 32

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 ) = R 3 V ( i ) ρ Euler ( t , x ; V ) dV , � � U ( i ) U ( j ) � ( t , x ) = R 3 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, or equivalently with a Lagrangian description of the flow. 33

Lagrangian description for turbulent flow On a given probability space (Ω , F , P ) , consider the fluid particle state vector ( X t , U t , θ t ) satisfying d X t = U t dt , � � − 1 d U t = ̺ ∇ x � P � ( t , X t ) − G ( t , X t ) ( U t − � U � ( t , X t )) + F t dt � + C ( t , X t ) ε ( t , X t ) dW t , d θ t = D 1 ( t , X t , θ t ) dt + D 2 ( t , X t , θ t ) d � W t . ( W , � W ) is a ( 3 d × 1 d ) -Brownian motion. � � U ( i ) � ( t , x ) , � U ( i ) U ( j ) � ( t , x ) are computed as conditional expectations. � ε , C , G , D 1 , D 2 are determined by a chosen Reynolds-averaged Navier-Stokes closure. 34

Compute the Reynolds averages � U ( i ) � and � U ( i ) U ( j ) � Let ρ Lagrangian ( t ; x , V , ψ ) be the probability density function of ( X t , U t , θ t ) . Eulerian PDF versus Lagrangian PDF: (case of incompressible flow, constant mass of particles) ρ Lagrangian ( t ; x , V , ψ ) � ρ Euler ( t , x ; V , φ ) dVd φ = R 4 ρ Lagrangian ( t ; x , U , ψ ) dUd ψ dVd φ For any random variable F ( U t , θ t ) , with finite moment, � � F ( U , ϑ ) � ( t , x ) = R 4 F ( V , ψ ) ρ Euler ( t , x ; V , ψ ) dVd ψ � � � � X t = x = E F ( U t , θ t ) In particular, � � � U ( i ) � X t = x � U ( i ) � ( t , x ) = E t � � � � U ( i ) U ( j ) � ( t , x ) = E U ( i ) U ( j ) � X t = x t t 35

Fluid particle model family On a given probability space (Ω , F , P ) , consider the fluid particle state vector ( X t , U t , ψ t ) satisfying d X t = U t dt � � − 1 d U t = ̺ ∇ x � P � ( t , X t ) − G ( t , X t ) ( U t − � U � ( t , X t )) + F t dt � + C ( t , X t ) ε ( t , X t ) dW t , d θ t = D 1 ( t , X t , θ t ) dt + D 2 ( t , X t , θ t ) d � W t . ( W , � W ) is a ( 3 d × 1 d ) -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 , D 1 , D 2 imposed by the chosen turbulence closure � compute ∇ x � P � � introduce boundary conditions 36

PDF description of flow via the associated Fokker Planck Equation P (( X t , U t , ψ t ) ∈ dxdvd φ ) := ρ Lagr. ( t ; x , v , φ ) dxdvd φ . ∂ t ρ Lagr. + ( v · ∇ x ρ Lagr. ) = 1 ̺ ( ∇ x � P � ( t , x ) · ∇ v ρ Lagr. ) − ∇ v · ( G ( t , x )( � U � ( t , x ) − v ) ρ Lagr. ) + 1 2 C ( t , x ) ε ( t , x ) △ v ρ Lagr. − ∂ φ ( D 1 ( t , x , φ ) ρ Lagr. ) + 1 2 ∂ 2 φφ ( D 2 ( t , x , φ ) ρ Lagr. ) � Integrating w.r.t. dvd φ gives the equation for the conservation of mass, � ̺ ( t , x ) = ρ Lagr. ( t ; x , v , φ ) dvd φ �� � � � v ρ Lagr. dvd φ � ∂ t ρ Lagr. dvd φ + ∇ x · ρ Lagr. dvd φ = 0 ⇐ ⇒ ∂ t ̺ + ∇ x · ( ̺ � U � ) = 0 . ρ Lagr. dvd φ � Integrating w.r.t. vdvd φ gives the RANS equation. � Integrating w.r.t. vv t dvd φ gives the Reynolds stress closure, according to C , G , D 1 , D 2 . 37

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,  with U t = ( u ( i ) i ( t ) = u ( i ) − � u ( i ) t , i = 1 , 2 , 3 ) and u ′ dX t = U t dt , t �  t  �� �� � � du ( i ) u ( j ) C 0 ε ( t , X t ) dB ( i ) − � u ( j ) � = − ∂ x i � P � ( t , X t ) dt + G ij ( t , X t ) dt +  t t t  j tke δ ij + C 2 ∂ � u ( i ) � ε G ij = − C R C 0 ε = 2 3 ( C R ε + C 2 P − ε ) , ∂ x j 2 P = 1 2 ( P 11 + P 22 + P 33 ) , the turbulent production term is computed with � � k � ∂ � u ( i ) � k � ∂ � u ( j ) � � � u ′ i u ′ + � u ′ j u ′ P ij := − . ∂ x k ∂ x k k ε is closed with mixing length or turbulent viscosity parametrizations. 38

The Reynolds tensor models in SLM Corresponding Reynolds Stress Equation : 3 ∂ t � u ′ i u ′ � U � · ∇ x � u ′ i u ′ � ∂ x k � u ′ i u ′ j u ′ � � j � + j � + k � = k = 1 � 1 � ε 3 2 + 3 � � � � u ′ i u ′ k � ∂ x k � U ( i ) � + � u ′ j u ′ k � ∂ x k � U ( j ) � tke � u ′ i u ′ − − 2 4 C 0 j � + C 0 εδ ij k = 1 versus Exact Reynolds Stress Equation : 3 ∂ t � u ′ i u ′ � U � · ∇ x � u ′ i u ′ � ∂ x k � u ′ i u ′ j u ′ � � j � + j � + k � = k = 1 3 3 − 1 � � j ∂ x i p ′ + u ′ � � u ′ i u ′ k � ∂ x k � U ( i ) � + � u ′ j u ′ k � ∂ x k � U ( j ) � ̺ � u ′ i ∂ x j p ′ � + ν � � ∂ x k u ′ i ∂ x k u ′ − j � k = 1 k = 1 � 3 j � = 2 � ∂ x l u ′ i ∂ x l u ′ The dissipation tensor is reduced to the relation ν 3 ε E δ ij l = 1 The turbulent pressure terms are closed according to the Rotta’s second order closure model : � � ε � � 1 1 + 3 j � + 2 � p ′ ∂ x j u ′ i � + � p ′ ∂ x i u ′ tke � u ′ i u ′ j � = − 3 C R εδ ij , 2 C 0 ̺ � � 1 + 3 for 2 C 0 denoting the prescribed Rotta’s constant. 39

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,  with U t = ( u ( i ) i ( t ) = u ( i ) − � u ( i ) t , i = 1 , 2 , 3 ) and u ′ dX t = U t dt , t �  t  �� �� � � du ( i ) u ( j ) C 0 ε ( t , X t ) dB ( i ) − � u ( j ) � = − ∂ x i � P � ( t , X t ) dt + G ij ( t , X t ) dt +  t t t  j tke δ ij + C 2 ∂ � u ( i ) � ε G ij = − C R C 0 ε = 2 3 ( C R ε + C 2 P − ε ) , ∂ x j 2 P = 1 2 ( P 11 + P 22 + P 33 ) , the turbulent production term is computed with � � k � ∂ � u ( i ) � k � ∂ � u ( j ) � � � u ′ i u ′ + � u ′ j u ′ P ij := − . ∂ x k ∂ x k k ε is closed with mixing length or turbulent viscosity parametrizations. 40

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]  with U t = ( u ( i ) i ( t ) = u ( i ) − � u ( i ) t , i = 1 , 2 , 3 ) and u ′  dX t = U t dt , t �    t  � � � � du ( i ) u ( j ) C 0 ε ( t , X t ) dB ( i ) − � u ( j ) �  ( t , X t ) dt +  = − ∂ x i � P � ( t , X t ) dt + G ij   t t t j ε G ij = − γ ij − 1 C 0 ε = 2 3 ( γ ij ) � u ′ i u ′ tke δ ij j � 2 homogeneous γ ij = ( 1 − α tke ) γ wall ij + α tke γ ij L 2 ∇ 2 α − α = − 1 Elliptic blending coefficent α solves: tke ( L length scale), and ∂ � u ( i ) � 2 ( C R − 1 ) ε = − 1 homogeneous − γ tke δ ij + C 2 ij ∂ x j = − 7 . 5 ε − γ wall tken i n j ij 41 where n is the wall-normal unit vector.

Dissipation closure in the atmospheric boundary layer collaboration with P . Drobinski (IPSL) With notation U t = ( u t , v t , w t ) and also U t − � U � ( t , X t ) = ( u ′ t , v ′ t , w ′ t ) � Closure relation classically used in meteorology (see e.g. Cuxart et al. 2000) : tke 3 / 2 ( t , x , y , z ) ε ( t , x , y , z ) = C ε . ℓ m ( z ) � Relation based on turbulent viscosity parametrization (O’Brien 70, Heinz 97) � � u ′ w ′ � 2 + � v ′ w ′ � 2 ν turb = �� � 2 + � � 2 ∂ z � u � ∂ z � v � ε = tke 2 1 1 � 1 � ν param 2 + 3 3 4 C 0 turb ℓ m ( z ) = κ z 1 [ 0 , z c ] + κ z c 1 [ z c , 750 ] 42

Wall-boundary condition The velocity is reflected at level z mirror (Drebeen Pope 98, Minier Cell center Cells • { z = zc } Pozorski 99) U in Outward particle • w in = − w out Inward normal u in = u out − 2 � u ′ w ′ � { z = z mirror } � w ′ 2 � w out , U out • Mirror particle • v in = v out − 2 � v ′ w ′ � { z = z 0 } � w ′ 2 � w out . { z = 0 } Ground We compute � u ′ w ′ � and � v ′ w ′ � using the MESO-NH-style formula u 2 u 2 ∗ � u � ∗ � v � � u ′ w ′ � ( t , x , y ) = − � v ′ w ′ � ( t , x , y ) = − � � u � 2 + � v � 2 ( t , x , y ) , � � u � 2 + � v � 2 ( t , x , y ) � � u � 2 ( t , x c , y c , z c ) + � v � 2 ( t , x c , y c , z c ) with u ∗ ( t , x c , y c ) = κ log ( z c / z 0 ) (see e.g. [Schmid and Schumann 89, Lafore et al 98, Carlotti 02]) 43

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

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 N p particles ( X k , N p , U k , N p ) in t t D . Each cell C contains N pc particles by constant mass density constraint. K ( ., x ) = 1 ( ., C ( x )) . N p � � � 1 N p � � 1 U k , N p K ( X k , N p K ( X k , N p � F ( U ) � ( t , x ) ≃ F , x ) , x ) t t t N p N p k = 1 k = 1 N p � K ( X k , N p , x ) = N pc t k = 1 45

The interacting particles system For j = 1 , . . . , N p  d X j , N p = U j , N p  dt ,  t t   = − 1 ∂ � P �  d U ( i ) , j , N p ( t , X j , N p ) dt t t ̺ ∂ x i   + D U ( t , X j , N p ) dt + B U ( t , X j , N p ) dW ( i ) , j , N p    t t t + jump terms for boundary condition , i ∈ { 1 , 2 , 3 } . • The coefficients D U , B U depend on the particles approximations of � U � , � U ( i ) U ( j ) � , and their derivatives. ∂ � P � • − 1 ( t , X j , N p ) ensures that ∇ · � U � = 0 and maintains the mass density t ̺ ∂ x i constant. ∇ 2 � P � = − ∂ 2 � U ( i ) U ( j ) � ∂ x i ∂ x j 46

Algorithm Fractional step method : n ∆ t − → ( n + 1 )∆ t (Pope 85, Minier Pozorski 99) Step 1 : For n ∆ t ≤ t ≤ ( n + 1 )∆ t , ( X j , N p n ∆ t , U ( i ) , j , N p , j = 1 , . . . , N p ) , given n ∆ t  X j , N p U j , N p d � = � dt ,   t t   ∂ � P � U ( i ) , j , N p X j , N p d � ( t , � = − 1 ) dt t t ̺ ∂ x i  U ( t , X j , N p U ( t , X j , N p ) dW ( i ) , j , N p  + D � ) dt + B �   t t t + jump terms , i ∈ { 1 , 2 , 3 } X j , N p → X j , N p Step 2: Correction of the particles positions � ( n + 1 )∆ t − ( n + 1 )∆ t that maintains the (discrete) uniform distribution, by solving a (discrete) optimal transport problem. U j , N p → U j , N p Step 3: Correction of the particles velocity � ( n + 1 )∆ t − ( n + 1 )∆ t such that ∇ · � U � ( n + 1 )∆ t = 0 , as in Chorin-Temam projection method : � ∆ t ∇ · � � 1 △ p = U � ( n + 1 )∆ t , x ∈ D , ∂ p ∂ n D = 0 , on ∂ D U j , N p U j , N p ( n + 1 )∆ t − ∆ t ∇ p ( X j , N p ( n + 1 )∆ t = � ( n + 1 )∆ t ) and � U � ( n + 1 )∆ t = � � U � ( n + 1 )∆ t − ∆ t ∇ p . 47

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

Actuator disc methods in the Lagrangian setting d X t = U t dt � � − 1 d U t = ρ ∇ x � P � ( t , X t ) dt � � − G ( t , X t ) U t − � U � ( t , X t ) dt + C ( t , X t ) dW t + f ( t , X t , U t ) dt + f nacelle ( t , X t , U t ) dt + f mast ( t , X t , U t ) dt . f ( t , X t , U t ) the body forces that the blades exert on the flow. Supplementary terms f nacelle ( t , X t , U t ) and f mast ( t , X t , U t ) represent the impact of the mill nacelle and mast. (B., Espina, Morice, Paris, Rousseau 14) 49

Forces Non rotating, uniformly loaded actuator disc model: knowing the axial induction factor a � � f x ( t , X t , U t ) = − 1 2 a � 2 e x . � E [ u t | X t ∈ C ] 1 − a ∆ x Rotating actuator disc : N b 4 π r ∆ xU relat ( X t , U t ) 2 c ( r ( X t )) ( C L ( α ) cos ( φ ) + C D ( α ) sin ( φ )) ( X t , U t ) , f x ( t , X t , U t ) = − 1 { X t ∈C} N b 4 π r ∆ xU relat ( X t , U t ) 2 c ( r ( X t )) ( C L ( α ) sin ( φ ) − C D ( α ) cos ( φ )) ( X t , U t ) . f θ ( t , X t , U t ) = 1 { X t ∈C} 50

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

Validation with wind tunnel measures Figure 2: Turbulence intensity with SDM, using the inflow mean velocity U hub at the hub � height given in Wu-Porte Agel 2011: I = 2 / 3 tke / U hub tke isosurface, for a (3 blades) actuator line model. 52

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

Wind farm evaluation 54

Wind farm evaluation 55

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

Well posedness of a simplified Langevin model Consider  � t  X t = X 0 + 0 U s ds , � t � t U t = U 0 + 0 B [ X s , U s ; ρ s ] ds + 0 σ ( s , X s , U s ) dW s ,  P (( X t , U t ) ∈ dxdu ) = ρ t ( x , u ) dxdu , t ∈ [ 0 , T ] , �  � R d b ( u , v ) γ ( x , v ) dv  � , if R d γ ( x , v ) dv � = 0 , where B [ x , u ; γ ] = R d γ ( x , v ) dv  0 , elsewhere If b : R 2 d → R d is bounded, by the Girsanov theorem, any weak solution ( X t , U t , t ∈ [ 0 , T ]) has a strictly positive density ( ρ t ( x , u ) , t ∈ [ 0 , T ]) , E [ b ( u , U t ) | X t = x ] = B [ x , u ; ρ t ] . and 57

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 ∈ R d , | v | 2 � d a ( i , j ) ( t , x , u ) v i v j ≤ λ | v | 2 . ≤ λ i , j = 1 ◮ For all 1 ≤ i , j ≤ d , a ( i , j ) is H¨ older continuous in the following sense: α ∈ ( 0 , 1 ] | a ( i , j ) ( t , x , u ) − a ( i , j ) ( s , y , v ) | α 2 + | x − y − v ( t − s ) | α + | u − v | α ) . ≤ C ( | t − s | Theorem Bossy, Jabir, Talay 11 � � Let b ∈ C b ( R 2 d , R d ) , let ( X 0 , U 0 ) s.t. E P � X 0 � R d + � U 0 � 2 < + ∞ . On the R d previous hypotheses on the velocity diffusion coefficient σ , the system has a unique weak solution. 58

The smoothed system in the space variables � t � X ε 0 U ε t = X 0 + s ds , � t � t U ε 0 B ε [ X ε s , U ε s ; ρ ε t = U 0 + s ] ds + 0 σ ( s , X s , U s ) dW s , where L aw ( X ε t , U ε t ) = ρ ε t ( x , u ) dxdu , and ∀ non–negative γ in L 1 ( R 2 d ) , � R 2 d b ( v , u ) φ ε ( x − y ) γ ( y , v ) dy dv � B ε [ x , u ; γ ] = R 2 d φ ε ( x − y ) γ ( y , v ) dy dv + ε , for a given regularization φ ε ∈ C 1 c ( R d ) of the Dirac mass in R d .  d X ε t = U ε t dt ,   �   �  E [ b ( v , U ε t ) φ ε ( x − X ε � t t )]  � x = X ε t , v = U ε d U ε σ ( s , X ε s , U ε � t = t dt + s ) dW s , (10) �   E [ φ ε ( x − X ε t )] + ε � 0    x = X ε  t t ∈ [ 0 , T ] . 59

Existence of the smoothed system On (Ω , F , ( F t ; t ∈ [ 0 , T ]) , Q ) , � � W i independent Brownian motions in R d , t ; t ∈ [ 0 , T ] ; N { ( X i 0 , U i 0 ); i ∈ N } i.i.d, independent of the Brownian family, such that � � ( X 1 0 , U 1 Q 0 ) ∈ dx du = ρ 0 ( x , u ) dx du .  � t X i ,ε, N 0 U i ,ε, N  = X i 0 + ds ,  t  s   � N   b ( U i ,ε, N , U j ,ε, N ) φ ε ( X i ,ε, N − X j ,ε, N  ) � t s s s s j = 1 , j � = i U i ,ε, N  ds + W i = U i � � 0 + t ,  � N  t  φ ε ( X i ,ε, N − X j ,ε, N  ) + ε 0  s s j = 1 , j � = i   i = 1 , . . . , N . � � � N The sequence { π N = L aw 1 i = 1 δ { X i ,ε, N , U i ,ε, N , W i } , N ∈ N } is tight on N P ( C ([ 0 , T ]; R 3 d )) . The sequence {L aw ( X i ,ε, N , U i ,ε, N , 1 ≤ i ≤ N ) , N } is L aw ( X ε , U ε ) –chaotic. 60

Convergence ε → 0 , uniqueness results The linear Langevin system (I) For ( y , v ) ∈ R 2 d , consider ( Y s , y , v , V s , v t ; t ≥ s ≥ 0 ) solution of the Langevin t equation � t   Y s , y , v V s , y , v  = y + d θ,  t θ s (11) � t   V s , y , v σ ( θ, Y s , y , v , V s , y , v  = v + ) dW θ . t θ s 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 ] × R 2 d , 1 ≤ i , j ≤ d , the derivatives ∂ v i Γ( s , y , v ; t , x , u ) exist and are continuous for all ( s , y , v ) ∈ R × R 2 d such that ( s , y , v ) � = ( t , x , u ) . (Ultraparabolic operators of Kolmogorov type (Di Francesco&Pascucci 05)) 61

Convergence ε → 0 , uniqueness results The linear Langevin system (II) Proposition (following) (ii) Let f ∈ C b ( R d ) . The function G t , f defined by G t , f ( s , y , v ) = E ( f ( Y s , y , v , V s , y , v )) , t t is the unique classical solution of the Cauchy problem:  ∂ s G t , f + A ( s , G t , f ) = 0 on [ 0 , t ) × R 2 d ,  s → t − G t , f ( s , y , v ) = f ( y , v ) ∀ ( y , v ) ∈ R 2 d , lim  � d A ( ψ )( s , y , v ) = ( u · ∇ x ψ ( s , y , v )) + 1 ( σσ ∗ ) ( i , j ) ( s , x , v ) ∂ 2 v i , v j ψ ( s , y , v ) 2 i , j = 1 (iii) There exists C > 0 depending only on T and λ such that � C R d |∇ v Γ( s , y , v ; t , x , u ) | dx du ≤ √ t − s , ∀ 0 ≤ s < t ≤ T . sup ( y , v ) ∈ R 2 d 62

Convergence ε → 0 , uniqueness results t , s ) is the adjoint of the transition semigroup of ( Y s , y , v , V s , y , v ( S ∗ ) , t t Define S ′ t , s : L 1 ( R 2 d ; R d ) → L 1 ( R 2 d ; R ) by � S ′ t , s ( h )( x , u ) = R 2 d ( ∇ 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 L 1 ( R 2 d ) � t S ∗ S ′ ρ t = t , 0 ( µ 0 ) + t , s ( ρ s ( · ) B [ · ; ρ s ]) ds 0 � t ρ ǫ S ∗ S ′ t , s ( ρ ǫ s ( · ) B ǫ [ · ; ρ ǫ = t , 0 ( µ 0 ) + s ]) ds t 0 For all t ∈ ( 0 , T ] , h , δ ∈ R d , it holds that � R 2 d | ρ ǫ t ( x + h , u + δ ) − ρ ǫ t ( x , u ) | dx du = 0 . | h | , | δ |→ 0 lim sup lim ǫ → 0 + 63

Definition - : The non linear Martingale Problem ( MP ε ) : A probability P ε on C ([ 0 , T ]; R 2 d ) is a weak solution of (10), or equivalently, a solution to the martingale problem ( MP ε ) if ( i ) P ε ◦ ( x 0 , u 0 ) − 1 = µ 0 . ( ii ) For all t ∈ ( 0 , T ] , the time marginal P ε ◦ ( x t , u t ) − 1 has a density ρ ǫ t w.r.t. Lebesgue measure on R 2 d . ( iii ) For all f ∈ C 2 b ( R 2 d ) , the process � t A ǫ f ( x t , u t ) − f ( x 0 , u 0 ) − s f ( s , x s , u s ) ds (12) ρ ǫ 0 is a P ε –martingale • Any limit point π ∞ of { π N , N ∈ N } has a full measure on the set of solutions of ( MP ε ) . • Consequently: Any subsequence of { π N , N ∈ N } converges to δ P ε or � � X i ,ε, N , U i ,ε, N , W i , i = 1 , . . . , N equivalently L aw is P ε -chaotic. 64

The limit ε → 0 Lemma : The limit Martingale Problem ( MP ) There exists P on C ([ 0 , T ]; R 2 d ) such that P is a solution to the martingale problem ( MP ) : ( i ) P ◦ ( x 0 , u 0 ) − 1 = µ 0 . ( ii ) For all t ∈ ( 0 , T ] , the time marginal P ◦ ( x t , u t ) − 1 has a positive density ρ t w.r.t. Lebesgue measure on R 2 d . ( iii ) For all f ∈ C 2 b ( R 2 d ) , the process � t f ( x t , u t ) − f ( x 0 , u 0 ) − A ρ s f ( s , x s , u s ) ds is a ( MP ) –martingale. 0 Main difficulty : for any bounded F s -measurable r.v Φ , � � � t ( B ε [ x θ , u θ ; ρ ε Φ θ ] · ∇ u f ( x θ , u θ )) d θ ε → 0 + E P ε lim s � � t � = E P Φ ( B [ x θ , u θ ; ρ θ ] · ∇ u f ( x θ , u θ )) d θ . s 65

Compute the wind at a finer scale Downscaling method for local wind assessment, based on a stochastic Lagrangian approach 66

The downscaling problem France 44 ◦ N 41 ◦ N Mer M´ editerran´ ee 0 ◦ 4 ◦ E 8 ◦ E Comparison WRF and SDM on flat terrain; horizontal visualization of the u component of the wind at the altitude 400 m ; time is T 0 + 4 days. 67

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 ) = V ext ( t , x ) where V ext is the WRF forcing. Bottom boundary condition (in blue) necessities to introduce a wall law model. 68

The boundary condition in the PDF approach ∀ x ∈ ∂ D , � U � ( t , x ) = V ext ( t , x ) . � R d v ρ ℓ ( t , x , v ) dv � R d ρ ℓ ( t , x , v ) dv = V ext ( t , x ) . � � R d v ρ ℓ ( t , x , v , ) dv = R d V ext ( t , x ) ρ ℓ ( t , x , v ) dv � � ⇔ R d v ρ ℓ ( t , x , v , ) dv = R d v ρ ℓ ( t , x , v + 2 ( V ext ( t , x ) − v )) dv ⇑ ρ ℓ ( t , x , v , ) = ρ ℓ ( t , x , v + 2 ( V ext ( t , x ) − v )) , ∀ v ∈ R d leads to specular boundary condition with jump on ∂ D for ρ ℓ ... 69

Stochastic Downscaling Method (SDM) Let D be an open set of R 3 , and a velocity V ext given at ∂ D : � t X t = X 0 + U s ds , 0 � t � � 1 � ε ( s , X s ) � − 1 2 + 3 d U t = U 0 + ̺ ∇� P � ( s , X s ) − tke ( s , X s ) ( U s − � U � ( s , X s )) 4 C 0 ds 0 � t � � + C 0 ε ( s , X s ) dW s + 2 ( V ext ( s , X s ) − U s − ) 1 { X s ∈ ∂ D} . 0 0 < s ≤ t The jump term should ensure that � U � ( t , x ) = V ext ( t , x ) , ∀ x ∈ ∂ D . Generic nonlinear SDE : � t X t = X 0 + 0 U s ds , � � t � t U t = U 0 + 0 B [ X s ; ρ ( s )] ds + 0 Σ[ X s ; ρ ( s )] dW s − 2 ( V ext − U s − ) 1 { X s ∈ ∂ D} , 0 < s ≤ t ρ ( t ) is the probability density of ( X t , U t ) for all t ∈ ( 0 , T ] , � � R d b ( v ) f ( x , u ) du R d σ ( u ) f ( x , u ) du � � � � B [ x ; f ] = R d f ( x , u ) du � = 0 } , Σ[ x ; f ] = R d f ( x , u ) du � = 0 } . 1 { 1 { R d f ( x , u ) du R d f ( x , u ) du 70

Well-posedness of the SLM with the no-permeability boundary condition ( � U � ( t , x ) · n D ( x )) = 0  � t X t = X 0 + 0 U s ds ,   � � t  U t = U 0 + 0 B [ X s ; ρ ( s )] ds + σ W t − 2 ( U s − · n D ( X s )) n D ( X s ) 1 { X s ∈ ∂ D} ,   0 < s ≤ t  ρ ( t ) is the probability density of ( X t , U t ) for all t ∈ ( 0 , T ] , Q T = ( 0 , T ) × D × R d Σ T = ( 0 , T ) × ∂ D × R d and 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 ∞ : Q t c � ( u · n D ( x )) γ ( ρ )( s , x , u ) f ( s , x , u ) ds d σ ∂ D ( x ) du Σ t ∂ s f + u · ∇ x f + B [ · ; ρ · ] · ∇ u f + σ 2 � � � � = D× R d ( f ( 0 ) ρ ( 0 ) − f ( t ) ρ ( t )) + 2 △ u f ρ ( s ) Q t and, for dt ⊗ d σ ∂ D a.e. ( t , x ) in ( 0 , T ) × ∂ D , � � R d | ( u · n D ( x )) | γ ( ρ )( t , x , u ) du < + ∞ , R d γ ( ρ )( t , x , u ) du > 0 . 71

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 ]; R d )) , and this law admits time marginal densities ρ t = P ◦ ( x ( t ) , u ( t )) − 1 in L 2 ( ω ; D × R d ) , for all t ∈ [ 0 , T ] , that solves in V 1 ( ω ; Q T )  ∂ t ρ ( t , x , u ) + u · ∇ x ρ ( t , x , u ) − σ 2   2 △ u ρ ( t , x , u ) = − ( B [ x ; ρ ( t )] · ∇ u ρ ( t , x , u )) in Q T   ρ ( 0 , x , u ) = ρ 0 ( x , u ) in D × R d ,     γ − ( ρ )( t , x , u ) = γ + ( ρ )( t , x , u − 2 ( u · n D ( x )) n D ( x )) in Σ − T , The trace γ ( ρ ) (in the sense of the previous definition) satisfies the no-permeability boundary condition � R d ( u · n D ( x )) γ ( ρ )( t , x , u ) du � E { ( U t · n D ( X t )) | X t = x } = = 0 , dt ⊗ d σ ∂ D R d γ ( ρ )( t , x , u ) du -a.e. on ( 0 , T ) × ∂ D . The associated smoothed N -particles system propagates chaos. 72

The hypotheses ( H ) ( H Langevin ) for the construction of the linear Langevin process ( H MVFP ) for the well-posedness of the Vlasov-Fokker-Planck equation ◮ ( H Langevin )-( i ) ( X 0 , U 0 ) is distributed according to the initial law µ 0 having its � � | x | 2 + | u | 2 � support in D × R d and such that µ 0 ( dx , du ) < + ∞ . D× R d ◮ ( H Langevin )-( ii ) The boundary ∂ D is a compact C 3 submanifold of R d . ◮ ( H MVFP )-( i ) b : R d → R d is a bounded continuous function. ◮ ( H MVFP )-( ii ) The initial law µ 0 has a density ρ 0 in the weighted space L 2 ( ω, D × R d ) with ω ( u ) := ( 1 + | u | 2 ) α for some α > d + 3 . 2 ◮ ( H MVFP )-( iii ) There exist two measurable functions P 0 , P 0 : R + − → R + such that 0 < P 0 ( | u | ) ≤ ρ 0 ( x , u ) ≤ P 0 ( | u | ) , a.e. on D × R d ; � 2 and R d ( 1 + | u | ) ω ( u ) P 0 ( | u | ) du < + ∞ . 73

Propagation of chaos With independent copies { ( X i 0 , U i 0 , ( W i )) , i = 1 , . . . , N } of ( X 0 , U 0 , ( W )) � t   X i ,ǫ, N = X i U i ,ǫ, N  0 + ds ,  t  s   0  � t  U i ,ǫ, N t + K i ,ǫ, N = U i B ǫ [ X i µ ǫ, N ] ds + σ W i 0 + s ; ¯ , t t s   0  � �  �   ) 1 � � , i = 1 , . . . , N . K i ,ǫ, N U i ,ǫ, N · n D ( X i ,ǫ, N n D ( X i ,ǫ, N  = − 2 )  X i ,ǫ, N t s − s s ∈ ∂ D s 0 < s ≤ t � � N D× R d b ( v ) β ǫ ( y ) φ ǫ ( x − y ) γ ( dy , dv ) µ ǫ, N = 1 � where ¯ δ ( X i ,ǫ, N ) , B ǫ [ x ; γ ] = , U i ,ǫ, N t N D× R d β ǫ ( y ) φ ǫ ( x − y ) γ ( dy , dv ) + ǫ t t i = 1 Theorem Let P be the law on E of ( X , U , K ) , and let P ǫ, N be the law of { ( X i ,ǫ, N , U i ,ǫ, N , K i ,ǫ, N ); 1 ≤ i ≤ N } . Then P ǫ, N is P -chaotic; namely, for all { F l , 1 ≤ l ≤ k } , k ≥ 2 , with F i ∈ C b ( C ([ 0 , T ]; D ) × D ([ 0 , T ]; R d ) × D ([ 0 , T ]; R d )) , � k N → + ∞ � F 1 ⊗ F 2 ⊗ · · · ⊗ F k ⊗ 1 ⊗ 1 ⊗ · · · , P ǫ, N � = lim lim � F l , P� . ǫ → 0 + l = 1 74

Incompressible stochastic Lagrangian model in the torus B., Fontbona, Jabin and Jabir 2013 Goal : well-posedness for the following process in T d × R d : � t X t = [ X 0 + U s ds ] modulo 1 0 � t � � β ( � U s � − α U s ) − ∇ P ( s , X s ) U t = U 0 + ds + σ W t . ρ ( s , X s ) 1 { ρ ( s , X s ) > 0 } 0 with L aw ( X t ) = ρ ( t , x ) dx , � U t � = E ( U t | X t ) and � � � d ∂ 2 E ( U i t U j ( t , x ) ∈ R + × T d . − △ x P ( t , x ) = t | X t = x ) ρ ( t , x ) , ij i , j = 1 75

Incompressible stochastic Lagrangian model in the torus Lemma Assume that the previous nonlinear SDE has a solution ( X , U ) such that � T � T � E | U t | 2 < ∞ ∀ t ∈ [ 0 , T ] , | U s | 2 ds < ∞ and E T d |∇ P ( s , x ) | dxds < ∞ . 0 0 Assume also that for all 1 − periodic function ϕ : R d → R of class C 1 we have E ( ∇ ϕ ( X 0 ) · U 0 ) = 0 . Then, a) E ( ∇ ϕ ( X t ) · U t ) = 0 for all t ∈ [ 0 , T ] , b) the process ( X t , t ∈ [ 0 , T ]) is stationary. Existence result for the Vlasov Fokker Planck equation ? 76

Local existence of analytical solutions  ∂ t f ( t , x , u ) + u · ∇ x f ( t , x , u ) = σ 2   2 △ u f ( t , x , u ) + β d f ( t , x , u ) + β u · ∇ u f ( t , x , u )     � � �    = 0 on ( 0 , T ] × T d × R d ,  + ∇ u f ( t , x , u ) · ∇ P ( t , x ) − βα R d vf ( t , x , v ) dv   f ( t = 0 ) = f 0 on T d × R d ,     �     R d v i v j f ( t , x , v ) dv , on [ 0 , T ] × T d , △ P ( t , x ) = − ∂ x i x j Theorem (d=1) Let ¯ λ > 0 and s ≥ d + 2 be an even integer. There exists κ 0 = κ 0 (¯ λ, s ) and λ, s ) ≥ 0 s.t. if f 0 : T d × R d → R of class C ∞ and T > 0 satisfy: r �→ κ 1 ( r , ¯ ◮ � � R d f 0 ( x , u ) du = 1 and ∇ x · R d uf 0 ( x , u ) du = 0 for all x ∈ T , s u f 0 � ∞ ≤ C 0 ( | k | + m )!( | l | + n )! ◮ � ( 1 + | u | 2 ) 2 D l x D k for some n , m ∈ N , all pair of | k | + | l | λ multiindexes k , l ∈ N N and a constant C 0 < κ 0 (¯ λ, s ) , and ◮ T < κ 1 ( C 0 , ¯ λ, s ) , then, a solution f of class C 1 , ∞ exists. 77

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

We simplify again the fluid particle model Goal : analyze the rate of convergence � t   X t = X 0 + U s ds   0 � t   U t = U 0 + B [ X s ; ρ s ] ds + W t , ρ t = L ( X t , U t ) , for all t ≤ T  0 � R d b ( v ) γ ( x , v ) dv R d × L 1 ( R d × R d ) ∋ ( x , γ ) �→ B [ x ; γ ] = � 1 { � R d γ ( x , v ) dv > 0 } R d γ ( x , v ) dv � ( t , x ) �→ B [ x ; ρ ( t )] = E [ b ( U t ) � X t = 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] � R d × R d b ( v ) K ǫ ( x − y ) ρ ( dy , dv ) B ǫ [ x ; ρ ] := , � R d × R d K ǫ ( x − y ) ρ ( dy , dv ) where K ǫ ( x ) := 1 ǫ d K ( x ǫ ) and K is a kernel function. � t  X ǫ U ǫ  t = X 0 + s ds   0 � t  U ǫ B ǫ [ X ǫ s , ρ ǫ ρ ǫ t = L ( X ǫ t , U ǫ  t = U 0 + s ] ds + W t , t ) , for all t ≤ T .  0 79

We simplify the fluid particle model Goal : analyze the rate of convergence � t  � �  X t = X 0 + U s ds mod 1   0 � t   U t = U 0 + B [ X s ; ρ s ] ds + W t , ρ t = L ( X t , U t ) , for all t ≤ T  0 � R d b ( v ) γ ( x , v ) dv T × L 1 ( T × R d ) ∋ ( x , γ ) �→ B [ x ; γ ] = � 1 { � R d γ ( x , v ) dv > 0 } R d γ ( x , v ) dv � ( t , x ) �→ B [ x ; ρ ( t )] = E [ b ( U t ) � X t = 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 � T × R d b ( v ) K ǫ ( x − y ) ρ ( dy , dv ) B ǫ [ x ; ρ ] := , � T × R d K ǫ ( x − y ) ρ ( dy , dv ) where K ǫ ( x ) := 1 ǫ d K ( x ǫ ) and K is a kernel function. � t  � � X ǫ U ǫ  t = X 0 + s ds mod 1   0 � t  U ǫ B ǫ [ X ǫ s , ρ ǫ ρ ǫ t = L ( X ǫ t , U ǫ  t = U 0 + s ] ds + W t , t ) , for all t ≤ T .  0 80

Numerical analysis of a non degenerate toy model (d=1) � t Y t = Y 0 + b [ Y s ; ρ s ] ds + σ w t , ρ t = L ( Y t ) . 0 Assume b symmetric b ( x , y ) = b ( x − y ) . b [ x , ρ t ] = E b ( x − Y t ) . Introduce N particles: Let ( w i , i = 1 , . . . , N ) N -Brownian motion, independent of the i.i.d. ( Y i 0 , i = 1 , . . . , N ) with law ρ 0 . We consider � t N N 1 � t = 1 � Y i , N = Y i b ( Y i , N − Y j , N s ) ds + σ w i µ N 0 + t , δ Y i , N . t s N N t 0 j = 1 i = 1 Lemma Assume that there exists r ≥ 2 such that b ( · ) is in C r + 1 ( R ) , and ρ 0 ∈ W r , ∞ ( R ) ∩ W r , 1 ( R ) . b Then ρ t ∈ W r , ∞ , and for all test function φ ∈ C 3 c ( R ) , there exists a positive constant C such that � � C � � φ, µ N � � E t − ρ t � � ≤ √ . sup N t ∈ [ 0 , T ] 81 where the constant C is of the form β exp ( α � b ′ � t ) .

Control of the drift error term N , ρ t ] − 1 � � � b [ Y 1 , N � b ( Y 1 , N − Y j , N ∀ t ∈ [ 0 , T ] , E drift ( t ) := E ) � . � � t t t N j = 1 � t Notice that, with Y i t = Y i 0 b [ Y i s ; ρ s ] ds + σ w i 0 + t ,   � t N s , ρ s ] − 1 � � t − Y 1 , N � − Y j , N E | Y 1  b [ Y 1 b ( Y 1 , N  ds | = E ) � � t s s � � N 0 j = 1 Since b is Lipschitz with constant L = � b ′ � ∞ , x �→ b [ x , ρ t ] is Lipschitz, and applying Gronwall Lemma � t � t � t t − Y 1 , N s − Y 1 , N E | Y 1 L E | Y 1 | ≤ | ds + E drift ( s ) ds ≤ exp ( L ( t − s )) E drift ( s ) ds . t s 0 0 0 N t , ρ t ] − 1 � � E drift ( t ) ≤ 3 L E | Y 1 , N � t − Y j − Y 1 � b [ Y 1 b ( Y 1 t | + E t ) � � t � N j = 1 � t N � E b ( x − Y t ) − 1 �� � � � b ( x − Y j ≤ 3 L exp ( L ( t − s )) E drift ( s ) ds + E t ) . � � � N x = Y 1 0 t j = 1 �� � 2 � � E b ( x − Y t ) − 1 ≤ C j = 1 b ( x − Y j The ( Y j , j � = 1 ) being i.i.d. E � � N � t ) N . � � N � x = Y 1 t Thus � t ≤ exp ( 4 Lt ) C exp ( − Lt ) C ≤ exp ( 3 Lt ) C E drift ( t ) ≤ 3 L exp ( L ( t − s )) E drift ( s ) ds + √ √ √ . N N N 0 82

A symptomatic case � t 1 Y t = Y 0 + 2 ρ s ( Y s ) ds + w t , ρ t = L ( Y t ) . 0 McKean 1967, Calderoni &Pulvirenti 1983, Sznitman & Varadhan 1986, Sznitman 1986. � t N � Y i , N = Y i K ǫ [ Y i , N s ; µ N s ] ds + σ w i µ N 0 + t , t = δ Y i , N t t 0 i = 1 � R K ǫ ( x − y ) γ ( dy ) with K ǫ ( z ) = 1 ǫ K ( z with K ǫ [ x ; γ ] := ǫ ) and K positive function of mass equal to 1. When ρ 0 ∈ W 2 , 1 ( R ) ∩ W 2 , ∞ ( R ) then ρ t ∈ W 2 , 1 ( R ) ∩ W 2 , ∞ ( R ) for all finite t > 0 and � ρ t − ρ ǫ t � L 1 ( R ) ≤ C ǫ 2 , � t where Y ǫ 2 K ǫ [ Y ǫ s ; ρ ǫ ρ ǫ t = L ( Y ǫ 1 t = Y 0 + s ] ds + w t , t ) . 0 � � � � � � φ, µ N ,ǫ − ρ t � ǫ 2 + exp ( C /ǫ ) � � From the previous lemma: sup t ∈ [ 0 , T ] E � ≤ C . √ N BURGERS : pas=0.1 / 16000 particules/ T= 0.2 /Diffusion= 1 / Epsilon = 0.500E-2 BURGERS : pas=0.1 / 16000 particules/ T= 0.2 /Diffusion= 1 / Epsilon = 0.100 0.7 0.6 solution exacte solution exacte Approximation Approximation 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0 0 Numerical experiments are more optimistic : -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 ( ρ 0 ( dx ) = 1 [ 0 , 1 ] with N = 16000 , T = 0 . 2 , σ = 1 , ǫ = 0 . 510 − 2 left, ǫ = 0 . 1 right) 83

Go back to McKean kinetic particle system ( X i t , U i t , t ≤ T , i = 1 , . . . , N ) are the strong solution of � t  � � X i X i U i t = 0 +  s ds mod 1   0 � t  U i t = U i B [ X i s , ρ s ] ds + W i ρ t = L ( X i t , U i  0 + t , t ) , for all t ≤ T ,  0 where we recall that ( X i 0 , U i 0 , i = 1 , . . . , N ) are i.i.d. with density ρ 0 . ( X i ,ǫ t , U i ,ǫ t , t ≤ T , i = 1 , . . . , N ) are the strong solution of � t  � � X i ,ǫ X i U i ,ǫ = 0 +  s ds mod 1  t  0 � t  U i ,ǫ B ǫ [ X i ,ǫ s , ρ ǫ ρ ǫ t = L ( X i ,ǫ t , U i ,ǫ = U i s ] ds + W i  0 + t , t ) , for all t ≤ T ,  t 0 � t  X i , N = X i U i , N  0 + s ds ,   t   0   � N  j = 1 b ( U j , N − X j , N � t 1 s ) K ǫ ( X i , N s ) s U i , N = U i N ds + W i 0 + t , for all t ≤ T . t � N j = 1 K ǫ ( X i , N − X j , N  1 s )   0 N    � �� �   Nadaraya-Watson estimator at X i , N s � � t ; ρ t ] − F ǫ [ X 1 , N � F [ X 1 ; µ N � Control the error seen by the particles : E t ] � . t 84

Rate of convergence For any f measurable bounded on R d , we set the kernel regression version of the conditional expectation : ( x , t ) �→ F [ x ; ρ t ] = E [ f ( U t ) | X t = x ] , � R d × R d f ( v ) K ǫ ( x − y ) ρ t ( dy , dv ) = E [ f ( U t ) K ǫ ( x − X t )] ( x , t ) �→ F ǫ [ x ; ρ t ] = � R d × R d K ǫ ( x − y ) ρ t ( dy , dv ) E [ K ǫ ( x − X t )] Assume that f and b smooth and bounded function with bounded derivatives; K > 0 Lipschitz bounded, with compact support. Theorem 1 For any 1 < p < 1 + 1 + 3 d and c > 0 , there exists a constant C such that for all α > 0 , 1 ǫ > 0 and N > 1 satisfying ≤ c , we bound the error seen by the particles by 1 α 2 ǫ d N p � � � � 1 1 1 1 t ; ρ t ] − F ǫ [ X 1 , N � F [ X 1 � ; µ N � t ] � ≤ C ǫ + αǫ d N + ǫ ( d + 1 ) p N + + . E 2 √ t 1 dp ǫ ( ǫ d N ) N p If we choose α = ǫ , the optimal rate of convergence is achieved for N = ǫ − ( d + 2 ) p and �� � � 1 − � F [ X 1 t ; ρ t ] − F ǫ [ X 1 , N ; µ N ( d + 2 ) p . E t ] ≤ N � t 85

Error decomposition � Smoothing error : there exists C , such that for all t ∈ [ 0 , T ] for all ǫ > 0 . E | F [ X t ; ρ 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 ,  2  � � � � � � � j � = i f ( U j ,ǫ t ) K ǫ ( z − X j ,ǫ � t ) 1 � �  � F ǫ [ z ; ρ ǫ  ρ ǫ ǫ d ( N − 1 ) + ǫ 2 t ] − t ( z ) dz ≤ C . E � � � j � = i K ǫ ( z − X j ,ǫ � t ) + ( N − 1 ) α T � 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 � N �� � t ) K ǫ ( X i , N j = 1 E ǫ 0 [ f ( U ǫ − X ǫ t )] � t C C � 0 , X j 0 , U j � � F ǫ [ X i , N ; ρ ǫ t ] − ≤ C ǫ + ǫ d + 1 N + αǫ d N . E � N � t 0 [ K ǫ ( X i , N j = 1 E ǫ − X ǫ t )] + N α 0 , X j 0 , U j t � The iteration error 86

The iteration error Step 1 : remove the denominators � N � N �� t ) K ǫ ( X i , N � j = 1 E ǫ 0 [ f ( U ǫ − X ǫ � t )] j = 1 f ( U j , N ) K ǫ ( X i , N − X j , N ) � 0 , X j 0 , U j t � t t t E � − � � N � N 0 [ K ǫ ( X i , N j = 1 K ǫ ( X i , N − X j , N j = 1 E ǫ − X ǫ t )] + N α ) + N α t t t 0 , X j 0 , U j �� � � N �� � ≤ C � t ) K ǫ ( X i , N t )] − f ( U j , N ) K ǫ ( X i , N − X j , N � E ǫ 0 [ f ( U ǫ − X ǫ N E ) 0 , X j 0 , U j t t t t j = 1 �� � N � � �� + C C 0 [ K ǫ ( X i , N t )] − K ǫ ( X i , N − X j , N � E ǫ − X ǫ � N E ) + 0 , X j 0 , U j t t t 1 ( ǫ d N ) p j = 1 Step 2 : control the derivative of u ( k ) t ,χ ( s , y , v ) = E s , y , v [ b ( k ) ( U ǫ t ) K ǫ ( χ − X ǫ t )] . 1 For all t ≤ T , and p such that 1 ≤ p < 1 + 1 + 3 d , we have � C T d χ |∇ v u T ,χ ( t , x , u ) | p p ≤ sup 2 d ( p − 1 ) . p 2 + 3 ( T − t ) x , u for all ǫ > 0 . In particular, this quantity is time integrable on [ 0 , T ] and the bound of the time integral is independent of ǫ . 87

The iteration error Step 3 : particle removal �� � � t � � � � in 1 � � ∇ v u ( k ) t ( s , X j , N s , U j , N B ǫ [ X j , N s ] − B ǫ [ X j , N s ; ρ ǫ s ; µ N s ) · s ] , N E � ds � t , X i , N 0 j � = i Let p > 1 be sufficiently small. There exists a constant C such that for all α , ǫ and N satisfying ǫ d N > 1 , � t � � � � |∇ v u ( k ) � B ( l ) � t ( s , X j , N s , U j , N ǫ [ X j , N s ; µ N s ] − B ( l ) ǫ [ X j , N s ; ρ ǫ E s ) | 1 · s ] ds t , X i , N 0 � t � � C | B ( l ) ǫ [ X j , N − 1 ; ρ ǫ s ] − B ( l ) ǫ [ X j , N − 1 ; µ N − 1 ≤ 2 E ] | ds p s s s ( t − s ) 0 � t � � N − 1 � � � C 1 � � K ǫ ( X k , N − 1 − X j , N − 1 ) − E ǫ 0 [ K ǫ ( X k , N − 1 − X ǫ + s )] E ds 0 , X j 0 , U j s s s 1 N − 1 α 2 ( ǫ d N ) p 0 j = 1 � � C 1 C + 1 + + αǫ d N , 1 1 ( ǫ d N ) α 2 ( ǫ d N ) p p for all t ≤ T and all j � = i . 88

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 ≤ c , 1 α 2 ( ǫ d N ) p we have N − 1 � � � � � K ǫ ( X k , N − 1 − X j , N − 1 ) − E ǫ 0 [ K ǫ ( X k , N − 1 − X ǫ 1 t )] E � � 0 , X j 0 , U j t t t N − 1 j = 1 N − 1 � � ≤ C � C � � K ǫ ( X k , N − X j , N s ) − E ǫ 0 [ K ǫ ( X k , N − X ǫ 1 s )] + sup E � � 0 , X j 0 , U j s s d N − 1 1 ǫ s ≤ t ( ǫ d N ) q p j = 1 � � | B ǫ [ X j , N − 1 ; ρ ǫ t ] − B ǫ [ X j , N − 1 ; µ N − 1 E ] | t t t � � | B ǫ [ X j , N ; ρ ǫ t ] − B ǫ [ X j , N ; µ N − 1 , N ≤ C E ] | t t t N − 1 � � � C C � � K ǫ ( X k , N − X j , N s ) − E ǫ 0 [ K ǫ ( X k , N − X ǫ 1 + C sup s )] + + αǫ d N , E � � s 0 , X j 0 , U j s N − 1 1 s ≤ t ( ǫ d N ) p j = 1 where µ N − 1 , N is the normalized empirical distribution of the N − 1 processes ( X k , N , U k , N ) for all k ≤ N − 1 : µ N − 1 , N = � N − 1 1 j = 1 δ { ( X j , N , U j , N ) } . N − 1 89

Smoothing error � � � ≤ E | F [ X t ; ρ t ] − F ǫ [ X t ; ρ ǫ � F [ X t ; ρ t ] − F ǫ [ X ǫ t ; ρ ǫ + E | F ǫ [ X t ; ρ ǫ t ] − F ǫ [ X ǫ t ; ρ ǫ E t ] t ] | t ] | . � �� � � �� � A B A : use the uniform density lower bound + kinetic semigroup properties B : use moreover that the kernel function K is in C 1 b ( R d ) 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 , |∇ F ǫ [ x ; ρ ǫ t ] | ≤ C . sup sup t ≤ T x ∈ T 90

Rate of convergence for the measured error For any f measurable bounded on R d , we set the kernel regression version of the drift � R d × R d f ( v ) K ǫ ( x − y ) ρ ( dy , dv ) F ǫ [ x ; ρ ] := , � R d × R d K ǫ ( x − y ) ρ ( dy , dv ) ( x , t ) �→ F [ x ; ρ t ] = E [ f ( U t ) | X t = x ] f and b smooth and bounded function with bounded derivatives Proposition 1 For any 1 < p < 1 + 1 + 3 d and c > 0 , there exists a constant C such that for all α > 0 , ǫ > 0 and N > 1 satisfying 1 ≤ c , 1 α 2 ǫ d N p we have � � � 1 1 1 1 � E | F ǫ [ X j , N ; ρ ǫ t ] − F ǫ [ X j , N ; µ N t ] ≤ C ǫ + αǫ d N + ǫ ( d + 1 ) p N + + . � 2 √ t t 1 dp ǫ ( ǫ d N ) N p If we choose α = ǫ , the optimal rate of convergence is achieved for N = ǫ − ( d + 2 ) p and �� � � 1 − � F ǫ [ X j , N ; ρ ǫ t ] − F ǫ [ X j , N ; µ N ( d + 2 ) p . E t ] ≤ N � t t 91

Regression Kernels We consider only boxed kernel K such that b 1 S 0 , r ( x ) ≤ K ( x ) ≤ 1 S 0 , R 92

Nonparametric regression estimate Let ( X i , U i ) be an i.i.d sequence of n random variables. Local averaging estimate : approximate F ( x ) = E [ f ( U ) | X = x ] with: n � F n ( x ) = W n , i ( x ) · f ( U i ) , i = 1 where the weight W n , i depends only on x and ( X j , f ( U j )) , j ≤ n . Kernel estimate : K ǫ ( x − X i ) � n W n , i ( x ) = j = 1 K ǫ ( x − X j ) , where K ǫ ( x ) = K ( x ǫ ) . Related to the minimization: � 1 K ǫ ( x − X i )( f ( U i ) − c ) 2 � � n F n ( x ) = arg min . n c i = 1 93

Nonparametric regression estimate Partitioning estimate : Let P n = { A n , 1 , A n , 2 , . . . } be a partition of the domain. Set 1 { X i ∈ A n , j } � n for x ∈ A n , j . W n , i ( x ) = , i = 1 1 { X i ∈ A n , j } Related to the least squares estimate: � 1 ( f ( U i ) − f n ( X i )) 2 � n � F n = arg min , n f n i = 1 where the arg min is taken on the set of piecewise constant functions adapted to P n . Also: local modeling estimates (e.g. local polynomial estimates), global modeling/least squares estimates, penalized modeling. . . 94

Convergence rate of the nonparametric estimates ◮ Consider the mean L 2 error: � � � � � � F ( x ) − F n ( x ) � 2 ρ ( x ) dx E . D ◮ Variance-Bias decomposition: � � � � � � � � 2 ρ ( x ) dx � F ( x ) − F n ( x ) | bias ( x ) | 2 ρ ( x ) dx . E = Var F n ( x ) ρ ( x ) dx + D ◮ The previous estimators are weakly universally consistent (Stone’s theorem). ◮ For ( p , C ) smooth conditional expectations, the lower optimal convergence 2 p rate is of order n − 2 p + d . Theorem For Lipschitz continuous regression functions, we have E � F n − F � 2 ≤ C 1 n ǫ d + dC ǫ 2 , for both partitioning estimate and kernel estimate (with naive kernel). 2 This gives an optimal rate of n − d + 2 . 95

Local averaging based particle algorithm Particle algorithm : 1. Initialization of ( X i 0 , U i 0 ) for 1 ≤ i ≤ n 2. Time loop : for t k = k ∆ t up to T : ◮ for each particle 1 ≤ i ≤ n : 2.1 update the position : X i , n = X i , n k − 1 + U i , n k − 1 k 2.2 add gaussian noise: U i , n = U i , n k − 1 + ( W i k − W i k − 1 ) k 2.3 drift computation : loop over each particle to calculate the W n , j ( X i , n k − 1 ) : j = 1 b ( U j , n k − 1 ) K n ( X i , n k − 1 − X j , n � N k − 1 ) U i , n k + = ∆ t j = 1 K n ( X i , n k − 1 − X j , n � N k − 1 ) ◮ End loop particles 3. Final mean field evaluation at x : Loop on all particles � n j = 1 f ( U j , n T ) K n ( x − X j , n T ) � n j = 1 K n ( x − X j , n T ) Complexity up to O ( n 2 ) for kernel estimates, and only O ( n ) for partitioning estimates. 96

Numerical simulation For d = 2 Test case dU t = Cdt − ∇ V ( X t ) dt + ( � U t � − 2 U t ) dt + dW t , � �� � � �� � “pressure” conditional drift with V ( x , y ) = 1 2 π cos ( 2 π x ) sin ( 2 π y ) − 1 2 x , X 0 = σ Z with L ( Z ) = N ( 0 , 1 ) , L ( U 0 ) = N ( 0 , 1 ) , σ = 0 . 3 . Reference simulation Regression estimate implemented with Epanechnikov kernel T = 2 , 128 times step, N = 10 5 particles, ǫ = ¯ 1 16 . 97

98

We want to observe the error of the scheme as � | F [ x ; ρ T ] − F ǫ [ x ; µ ǫ, N T ] | ρ T ( x ) dx , D First we approximate the mean field F [ x ; ρ T ] by N mc � � � � 1 x ; µ ǫ, ¯ x ; µ ǫ, ¯ N ( ω k ) N ] := , F ǫ F ǫ N mc k = 1 Using spline interpolation on a mesh of size ∆ x , we observe the approximate L 1 error � � N mc � � � 1 ∆ x � F ǫ x ; µ ǫ, ¯ � ¯ − F ∆ x ǫ [ x ; µ ǫ, N ( ω k )] ρ ∆ x ( x ) dx . N ] N mc D k = 1 99

Observe the approximate L 1 error: 0.028 10 5 3 0.016 0.021 Optimal N ( ǫ ) 0.015 1 0 0 N = 1 /ǫ 4 0.014 . . 0 0 0.029 3 0.024 1 0.023 0.018 0.015 0.017 N 0.025 0.020 0.022 0.019 0.026 0.026 0 0.022 0.024 0.023 . 0 2 5 10 − 1 ǫ in-line with what we expected from the theoretical rate where the optimal value 1 of window size is given by d + 2 . 1 N 100