The Moment Guided Monte Carlo Method

The Momen t Guided Mo n te Carlo Metho d ∗ Pierre Degond 1,2 , Giacomo Dimarco † 1,2,3 and Lorenzo P aresc hi 4 1 Univ ersit ´ e de T oulouse; UPS, INSA, UT1, UTM ; Institut de Math ´ ematiques de T oulouse ; F-31062 T oulouse , F rance. 2 CNRS; Institut de Math ´ ematiques de T oulouse UMR 5219; F-31062 T oulouse, F rance. 3 Commissariat ` a l’Energie A tomique CEA-Sacla y DM2S-SFME; 91191 Gif-sur-Yv ette, F rance. 4 Univ ersit y of F errara; Departmen t of Mathematics; 44100 F errara, Italy . Octob er 22, 2018 Abstract In this work w e pro po se a new approa ch for the numerical sim ula tio n of kinetic equations through Monte Carlo schemes. W e introduce a new technique which p ermits to reduce the v ariance of particle metho ds through a ma tching with a set of s uitable macrosco pic mo ment equations. In or der to guar antee that the moment equations provide the correct so lutions, they ar e coupled to the kinetic equation throug h a non equilibrium ter m. The basic idea, o n which the method relies, co nsists in guiding the particle positio ns and velocities through momen t equations so that the concurrent solution of the mo ment and kinetic models furnishes the same macrosco pic quantities. Keyw ords: Mon te Carlo metho ds, h yb rid metho ds, v ariance reduction, Boltzmann equation, fl uid equations. Con ten ts 1 In tro duction 2 ∗ Ackno wledgements: This work w as supp orted by the Marie Curie Actions of the Europ ean Commission in the frame of the DEASE pro ject (MEST-CT-2005-02112 2) a n d by the F renc h Commisariat ` a l’ ´ Energie Atomique (CEA) in the frame of the contract ASTR E ( S A V 34160) † Corresponding author a d d ress: In stitu t de Math´ ematiqu es de T oulouse, UMR 5219 Univers it´ e Paul Sabatier, 118, route de N arb onne 31062 TOULOUSE Cedex, FRA NCE. E-mail addr esses :pierre.dego nd @math.univ-tou louse.fr, giacomo .dimarco@unife.it,lorenzo.paresc hi@unife.it 1 2 The Boltzmann equat ion and its fluid limit 3 3 The Momen t Guided Mon te Carlo Metho ds 4 3.1 Solution of the Momen t Equations . . . . . . . . . . . . . . . . . . . . . . . 5 3.2 The Momen t Matc hing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 4 Numerical results 9 5 Conclusions 14 1 In tro du ction The Boltzmann equation pro vides a kinetic description of gase s and more generally of particle systems. In many applications, the correct physical solution for a system far from thermo dyn amical equilibrium , such as, for instance, r arefi ed gases or plasmas requires the resolution o f the Boltzmann equati on [6]. The n umerical simulat ion of th e Bol tzmann equation with d eterministic tec hniqu es presents several drawbac ks due to difficulties i n treating the collision terms and to the large dimension of the pr oblem. The distribution function d ep ends on seven indep enden t v ariables: three co ordinates in physical space, thr ee co ordinates in v elo city s pace and the time. As a consequence, probabilistic tec hn iques s uc h as Direct Simulatio n Mon te Carlo (DSMC) metho d s are extensiv ely used in r eal situations due to their large fl exibilit y an d lo w computational cost compared to finite volume, fin ite difference or sp ectral methods for kinetic equ ations [1, 4, 5, 18]. On the other hand DSMC solutions are affecte d by large fl uctuations. Moreo v er, in n on stationary situations it is imp ossible to use time av erage s to reduce fluctuations and this leads to, either p o orly accurate solutions, or compu tationally exp ensiv e simulatio ns. More generally Mon te Carlo metho d s are frequentl y used in many real app lications to sim ulate ph ysical, chemica l and mathemat ical systems [16]. W e quote [4] for an o v erview on efficien t and low v ariance Mon te Carlo m etho ds. F or applications of v ariance reduction tec hniques to kinetic equation w e m en tion the works of Homolle and Hadjiconstantinou [12] and [13]. W e men tion also th e w ork of Boyd and Burt [3] and of P ullin [22] whic h dev elop ed a lo w diffusion particle metho d for sim ulating compressible in viscid flo ws. W e finally quote the w orks of Dimarco and P aresc hi [9] whic h w orked on the construction of efficien t and lo w v ariance metho ds for k in etic equations in transitional regimes. The basic idea describ ed in this work consists in reducing the v ariance of Monte C arlo metho ds by forcing particles to matc h prescrib ed sets of moment s giv en by the solution of deterministic equations. In order to pro vide th e correct solution, the moment equations are coupled to th e DSMC simulatio n of the Boltzmann equatio n through a kin etic correction term, wh ic h take s in to ac count dep artu res from thermo dynamical equilibrium. W e remark that the general method ology describ ed here is in dep end en t f r om the c hoice of the collisional kernel (Boltzmann, F okk er-Planc k, BGK etc..). Ho w ev er w e p oin t out that additional improv ements can b e obtained with hypotheses on the structure of the dis- tribution function, on the t yp e of considered k ernel and on the t yp e of resolution metho ds used for the kinetic and fluid equations. 2 In the p resen t pap er w e will f o cus on th e basic matc hin g tec hn ique, wh ic h consists in matc h ing the kinetic solution to that obtained by the deterministic solution of the first three momen t equations. The idea is that th e deterministic solution of the momen t equation (through finite vo lume or fin ite difference tec h niques) leads to a more accurate so- lution, in term of statistic al fl uctuations, than the DSMC method. Th erefore, w e co n s train the DSMC metho d to matc h the momen ts obtained through the deterministic resolution of th e m oment equations in such a w a y that th e higher accuracy of the m omen t resolution impro ves the accuracy of the DSMC metho d. W e exp erimen tally sho w that this is indeed the case. W e lea ve an in depth discussion of p ossible higher order matc hing extensions to futur e w ork. F or simplicit y , in the numerica l tests, we will make u se of a BGK collision term. Ho w eve r the form ulation of the m etho d is general and extensions to the full Boltzmann in teraction term are p ossib le without changing the stru cture of the algorithm as explained in details in the p ap er. Results in the case of the Bol tzmann op erator and impro v ements of the b asic tec hn iqu e d escrib ed here will b e present ed in [7 ]. The remaind er of the pap er is organized as f ollo ws. In th e next section we recall some basic notions on the Boltzmann equations and its fl uid limit. T he details of the numerical metho d are describ ed in Section 3. In Section 4 numerical examples w hic h demonstrate the capabilit y of the method are pr esen ted. Finally some futu re deve lopments and remarks are detailed in the last Section. 2 The Boltzmann equation and its flu id limit W e consider equ ations of the follo wing form ∂ t f + v · ∇ x f = Q ( f , f ) (1) with initial data f | t =0 = f init (2) where f ( x, v, t ) is a non negativ e function d escribing the time ev olution of th e distribu tion of particle s with v elo cit y v ∈ R d v and position x ∈ Ω ⊂ R d x at time t > 0. The operator Q ( f , f ) describ es particles inte r actions and is assumed to satisfy the lo cal conserv ation prop erties h mQ ( f , f ) i = 0 (3) where we define in tegrals o v er the ve lo city space as follo w s Z R d v ψ dv =: h ψ i (4) and m ( v ) = (1 , v , | v | 2 2 ) are the collision inv ariant s. Integrating (1) against its inv arian ts in v elo cit y space leads to the follo wing set of co nser v ations laws ∂ t h mf i + ∇ x h v mf i = 0 . (5) 3 Equilibrium fu nctions for the op erator Q ( f , f ) (i.e. s olutions of Q ( f , f ) = 0) are lo cal Maxw ellian of the form M f ( ρ, u, T ) = ρ (2 π T ) d v / 2 exp  −| u − v | 2 2 T  , (6) where ρ , u , T are the densit y , mean v elo cit y and temp erature of the gas at p osition x and at time t ρ = Z R d f dv , u = 1  Z R d v f dv , T = 1 d Z R d | v − u | 2 f dv . (7) In the s equ el we will den ote b y U = ( ρ, u, T ) , E [ U ] = M f . (8) Clearly we ha ve U = h mE [ U ] i . (9) No w, w hen the mean free path b et w een th e particles is very small compared to the t ypical length scale of the exp er im ent, the op erator Q ( f , f ) is large and we can r escale the space and time v ariables in (1) as x ′ = εx, t ′ = εt, (10) to obtain ∂ t f + v · ∇ x f = 1 ε Q ( f , f ) , (11) where ε is a small parameter p r op ortional to th e mean fr ee path and the p rimes ha ve b een omitted to k eep notations simple. P assing to the limit f or ε → 0 leads to f → E [ U ] and thus w e h a v e a closed h yp erb olic system of equat ions for the m acroscopic v ariables U ∂ t U + ∇ x F ( U ) = 0 , (12) with F ( U ) = h v mE [ U ] i . 3 The Momen t Gu ided Mon te Carlo Metho ds F or the sake of sim p licit y , in this work, we consider th e problem in one d im en sion b oth in physic al and v elo cit y spaces. Extensions to m ultidimensional problems are straigh tforw ard and will b e considered in [7]. The s tarting p oint of the method is the follo win g micro-macro decomp osition f = E [ U ] + g . ( 13) The function g represen ts th e n on-equilibrium part of the distribu tion f unction. F rom the definition ab ov e, it follo ws that g is in general n on p ositiv e. Moreo v er since f and E [ U ] ha ve the same momen ts w e ha ve h mg i = 0 . (14) 4 No w U and g satisfy the coupled system of equations ∂ t U + ∂ x F ( U ) + ∂ x h v mg i = 0 , (15) ∂ t f + v ∂ x f = Q ( f , f ) . (16) W e skip th e elemen tary p ro of of the ab o ve statemen t and r efer to [8] for details on the decomp osition of the distribution fu nction and the coupled systems whic h it is p ossible to deriv e. Our goal is to solv e the kinetic equation with a Mon te Carlo metho d, and concurr en tly the fluid equation with any type of finite difference or finite volume sc heme, where the correction term ∂ x h v mg i is ev aluated using particle moments. The tw o equations (15 -16), except for numerical errors, giv e the same r esults in terms of macroscopic quanti ties. It is natural to assume that th e set of m omen ts obtai n ed from the fluid sy s tem represents a b etter statistical estima te of the tru e m omen ts of the solution, since the resolution of the momen t equatio ns do es not in v olv e an y stoc hastic pro cess. Th u s we can summarize the method in the follo wing w a y . A t eac h time step t n 1. Solv e the kinetic equation (16) w ith a Monte Carlo sc h eme and obtai n a fir s t set of momen ts U ∗ = h mf ∗ i . 2. Solv e the fluid equation (15) with a fin ite v olume/difference scheme using particles to ev aluate ∂ x h v mg i and obtain a seco n d set of m omen ts U n +1 . 3. Matc h the momen ts of the kineti c solution with the fluid solution through a trans - formation of the samp les v alues f n +1 = T ( f ∗ ) so that h mf n +1 i = U n +1 . 4. Restart the c ompu tation to the n ext time step. F or Step 1, one can use any Monte Carlo metho d (or m ore generally an y lo w accurate but fast solve r ). Step 2 and 3 of the ab ov e procedu re requir e great care since they in vol ve the ev aluation of ∂ x h v mg i and the moment matc hing procedu re. Finally let us note that, in prin ciple, it is p ossible to impro v e the metho d, adding to system (15) additional equations for the time ev olution of higher ord er momen ts and get ∂ t h m n f i + ∂ x h v m n f i = h m n Q ( f , f ) i (17) with m n = v n and n ≥ 3. The solution of (17) with a finite vo lum e/difference sc heme, whic h in the general case is not straigh tforwa r d , will provide a b etter estimate of th e momen ts whic h are used in the momen t matc hin g [10], [14], [19 ]. W e will call this general class of metho ds Moment Guided Monte Carlo sc hemes. In the sequel, w e briefly focus o n steps 2 and 3 of the ab o ve p ro cedure. 3.1 Solution of the Momen t Equations In this section we discuss th e discretization of the momen t equations. W e will, at the end of the section, su ggest some appr oac hes that can p ossib ly b e u sed to impr o v e th e metho d 5 in the nearby future. Our scop e, in the construction of the numerical sc h eme, is to tak e adv antag e from the kno wledge of the Euler part of the momen t equations ∂ t U + ∂ x F ( U ) | {z } Euler equations + ∂ x h v mg i = 0 . (18) Th u s, the metho d is based on solving first the set of compressible Euler equations, for whic h many efficient n umerical metho ds ha ve b een dev elop ed in the literature, and then considering th e discretization of the kin etic fl ux ∂ x h v mg i . T o that aim, for the sp ace d is- cretizatio n of the compressible Euler equations we u s e b oth a first ord er central sc heme of Lax-F riedr ic hs t yp e or a second ord er MUSCL cent ral sc heme, w h ile a backw ard d is- cretizatio n is used for the time d eriv ative in all c ases U ∗ i − U n i ∆ t + ψ i +1 / 2 ( U n ) − ψ i − 1 / 2 ( U n ) ∆ x = 0 . (19) In the case of the second order s c heme the discrete flu x reads ψ i +1 / 2 ( U n ) = 1 2 ( F ( U n i ) + F ( U n i +1 )) − 1 2 α ( U n i +1 − U n i ) + 1 4 ( σ n, + i − σ n, − i +1 ) (20) where σ n, ± i =  F ( U n i +1 ) ± αU n i +1 − F ( U n i ) ∓ αU n i  ϕ ( χ n, ± i ) (21) with ϕ a give n slop e limiter, α equal to the larger eige nv alue of the Euler sys tem and χ n, ± i = F ( U n i ) ± αU n i − F ( U n i − 1 ) ∓ αU n i − 1 F ( U n i +1 ) ± αU n i +1 − F ( U n i ) ∓ αU n i (22) where the ab o v e ratio of v ectors is defined comp onen t wise. In the numerica l test section w e used b oth first and second order d iscretizatio n tec hn iques to test their differen t b eha viors when coupled with the DSMC solv er. As explained in that section, second order sc hemes increase fluctuations while fir st ord er do es not. Moreo v er, sin ce the numerical diffusion in tro d uced by first order schemes can b e excessiv e, we also test a switc hing metho d wh ich passes fr om first to second order accuracy as a fun ction of the r atio b et w een the n on equilibrium term and the equilibrium one. Thus, w e defi n e the foll owing quan tities β n i = λ n i | F 3 ( U n i ) | , λ i = Z R 3 v | v | 2 2 g n i dv (23) where F 3 ( U ) is the energy flu x . Note that in the one dimensional case the non equ ilibrium mass and momentum fluxes are identic ally zero. The MUSCL second ord er sc heme is then used when β n i is small while th e first order sc heme is used otherwise. This metho d, as sho wed in the tests, do es not increase flu ctuations and, at the same time, guarantee s a lo w er lev el of numerical dissipation in the results. W e n o w discuss ho w t o d iscretize the non equilibrium term ∂ x < v mg > . T o this aim, the same space first ord er discrete deriv ativ e is used as for the hydrod y n amic flu x F ( U ). The non equilibrium term h vmg i = h v m ( f − E [ U ]) i is c ompu ted by taking the difference 6 b et ween the momen ts of the particle solution and those of the Maxw ellian equilibrium. Th u s the fin al sc heme, for the moment equations, reads U n +1 i − U n i ∆ t + ψ i +1 / 2 ( U n ) − ψ i − 1 / 2 ( U n ) ∆ x + Ψ i +1 / 2 ( < v mg n > ) − Ψ i − 1 / 2 ( < v mg n > ) ∆ x = 0 . (24) where ψ i +1 / 2 ( U n ) can b e either first or second order while Ψ i +1 / 2 ( < v mg n > ) is alwa ys first order. In our metho d , in addition to the m ass momentum and energy equations, w e consider the third order moment ev olution equ ation. In this case, like in (17), we ha ve the additional problem of ev aluating the source term that now app ears at the righ t hand side. A t the particle lev el this can b e done b y simply measuring the v ariations of higher order momen ts in eac h cell during p article co llisions. Thus the discretized third o r d er momen t equation is p erformed in tw o steps, where the second o n e reads < m 3 f n +1 i > − < m 3 f ∗ i > ∆ t + Ψ i +1 / 2 ( < v m 3 f n > ) − Ψ i − 1 / 2 ( < v m 3 f n > ) ∆ x = 0 , (25) and w h ere f ∗ is the solution of the first step, th e collision step, wh ic h dep ends on the t yp e of collisional op erator. W e w ill describ e this step in the n umerical test section b elo w in the case of a BGK t yp e k ern el. The m ain adv an tage of considering additional moment equa- tions is that this red u ces the flu ctuations in the ev aluations of th e macroscopic quan tities U . Indeed, in the extended momen t system, particles pla y a role only in the ev aluation of higher order terms h v p f i , p > 3 and not dir ectly on the ev olution of the hydro dynamics quan tities. As a conclusion for this section, we d iscuss some p ossible improv ements whic h will b e dev elop ed in fu ture w orks [7 ]. T o this aim, w e observ e th at the d ecomp osition of the fl u x term int o an equilibrium and a non equilibrium part can be further exploited. Ind eed, as an effect of the guided Mon te Carlo tec hnique, the only remainin g sour ce of fluctuations in the moment equations is du e to the non equilibrium term ∂ x h v mg i . Thus, instead of u sing the same numerical sc heme as for the flu x ∂ x F ( U ), we can dev elop a sp ecific discretization metho d wh ic h further reduces the v ariance of these fluctuations. W e can consid er cell a v erages of the f orm 1 ∆ x Z x i +1 / 2 x i − 1 / 2 ∂ x h v mg i dx = h v mg i| x = x i +1 / 2 − h v mg i| x = x i − 1 / 2 ∆ x , m = (1 , v , | v | 2 ) , (26) with h v mg i = h v mf i − F ( U ). The int egral ov er the v elo cit y space can b e ev aluated b y summing ov er the particles h v mf i| x = x i +1 / 2 ≈ 1 N X j ∈ I i +1 / 2 B ( p j − x i +1 / 2 ) m j (27) where p j and ν j represent the p osition and velocit y of the j-th particle, m j = (1 , ν j , | ν j | 2 ), I i +1 / 2 a giv en space in terv al of size h (t ypically h ≥ ∆ x ) con taining x i +1 / 2 and B ≥ 0 is a suitable w eigh t function s.t. Z R B ( x ) dx = 1 . 7 F or example B ( x ) = 1 /h if | x | ≤ h/ 2 and B ( x ) = 0 elsewhere, giv es rise to a simp le sum of th e particles moments in the interv al I i +1 / 2 kno wn as the ’Nearest Grid Poin t’ pro cedur e in plasma physics [2]. S mo other reconstructions can b e reco vered by co nv olving the samples with a b ell-shap ed weig ht lik e a B-spline [21]. Note th at the v alue h has a strong influ ence on the fluctuations in the reconstructed function, an d in general should b e selected as a goo d compromise b etw een fluctuations and resolution. 3.2 The Moment Matc hing In the present work w e restrict ourselve s to the follo wing linear transformation: let a set of v elo cities ν 1 , . . . , ν J with fir st t w o moment s µ 1 and µ 2 b e giv en. Su pp ose b etter estimates σ 1 and σ 2 of the same moments are av ailable (usin g th e m oment equation). W e can apply the transf ormation describ ed in [4] ν ∗ j = ( ν j − µ 1 ) /c + σ 1 c = s µ 2 − µ 2 1 σ 2 − σ 2 1 , i = 1 , . . . , J (28) to get 1 J J X j =1 ν ∗ j = σ 1 , 1 J J X j =1 ( ν ∗ j ) 2 = σ 2 . Of course this renormalization is not p ossible for the moment of order zero (the mass densit y). Let us denote b y µ 0 an estimate of the zero order momen t and b y σ 0 its b etter ev aluation b y the momen t equations. Among the p ossib le techniques that can b e used to r estore a prescrib ed density we c ho ose to r eplicate or discard p articles inside the cells. Other p ossibilities are to deal with w eigh ted particles, mo v e particles among c ells according to some in terp olati on pro cedu r e or reconstruct the pr obabilit y distrib ution starting from samples and resample particles. W e lea ve a d eep er analysis of p ossible alternate choic es to futu re wo rk s . In order to reco v er the moment σ 0 , in the case µ 0 > σ 0 , we can use a discarding pro cedur e. Note that we w ould lik e to eliminate exactly th e f ollo wing num b er of p articles, f N p = µ 0 − σ 0 M p (29) where M p is the mass of a single p article. In general, the precise matc h is imp ossible, since the p articles mass is k ept fixed in time, and f N p can n ever b e an intege r. A fixed mass M p implies that µ 0 = N 1 M p , σ 0 = N 2 M p (30) with N 1 and N 2 in tegers suc h that N 1 > N 2 . N 1 and N 2 are the num b er o f particles in the cell b efore and a fter the matc h in g. Moreo v er, since the estimate σ 0 is not in ge neral an int eger m ultiple of M p , a mismatc h e such that e < ± M p is unav oidable. Thus w e can simply eliminate from the cell a su itable sto chasti c in teger appro ximation of f N p N p = Iround  µ 0 − σ 0 M p  (31) 8 where Iroun d( x ) is a stochastic roundin g d efi ned as Iround =  ⌊ x ⌋ + 1 , with probabilit y x − ⌊ x ⌋ ⌊ x ⌋ , with pr ob ab ility 1 − x + ⌊ x ⌋ (32) with ⌊ x ⌋ the intege r part of x . In the opp osite case, in which the mass of the particles insid e a cell is lo wer than the mass pr escrib ed b y the fluid equ ations µ 0 < σ 0 , the situation is less simple. In this situation, since the distr ib ution fu nction is not known analytically , it is not p ossible to sample n ew p articles w ithout introd ucing corr elations b et wee n samples. In this case we need to r eplicate N p = Iround  σ 0 − µ 0 M p  (33) particles. Note that this is d one all owing rep etitions. After the generati on step, samples are relo cated uniformly inside eac h spatial cell. No w we briefly discuss the p ossibilit y of f orcing samples to follo w higher ord er pre- scrib ed momen ts. T o this aim, obser ve that, the moment matc hing pr o cedure h as in fi nite p ossible solutions, since the num b er of particles inside a cell is larger than the num b er of the constraint s. Ho wev er, we aim at findin g a transformation w hic h p ossibly preserve s the Gaussian d itribution. T he only op erations whic h ob ey this constraint are linea r transfor- mations like (28), i.e. shifts a nd homotheties of the particle v elo cities. Ho w eve r, if w e slightl y relax the constrain t of preserv ation of the Gaussian d istr ibution, w e can reform u late the problem in the follo wing terms : find a suitable transform ation whic h leads to th e requ ir ed momen ts with the min imal change s in the d istribution f unction. In the general case, this request has a non trivial answ er whic h can b e reco vered b y solving an appropriate non linear system of e qu ations with sev eral constrain ts at eac h time step for eve ry cell. F or this reason an efficien t implementat ion of this pr o cedure is still an op en question. 4 Numerical results In the pr esent section w e r ep ort on some numerical results of th e moment guided metho d on different test cases obtained u sing a simplified BGK mo del for th e kinetic equation. First, we p erf orm an accuracy test usin g a smo oth p erio d ic solution and then we consider t w o classica l sho ck problems. In all the tests, we compare the momen t gu id ed (MG) solution with the standard Mon te Carlo (M C) solution and with the direct deterministic solution to the BGK equations based on a discrete vel o cit y mo del (DVM ) [17]. The Moment Guided DSMC metho d applied to t he BGK mo del In this paragraph we detail a p ossible algorithm, whic h m erges the tec hniqu es describ ed in the previous sections, in th e case of the simplified BGK collision operator. As usual the starting p oint of Mon te Carlo metho ds is giv en by a time sp litting [21] b et ween free tr ansp ort ∂ t f + v · ∇ x f = 0 , (34) 9 and collision, whic h in the case of the BGK op erator is substituted by a relaxation to wa rd s the equilibr ium ∂ t f = 1 ε ( f − E [ U ]) . (35) In Mon te Carlo sim ulations the distribu tion function f is discretized by a finite set of particles f = N X i =1 M P δ ( x − x i ( t )) δ ( v − v i ( t )) , (36) where x i ( t ) represent s th e particle p osition and v i ( t ) the particle velocit y . During the transp ort ste p then, the particles mo ve to their next p ositio ns according to x i ( t + ∆ t ) = x i ( t ) + v i ( t )∆ t (37) where ∆ t is such that an appropriate CFL condition holds. The collision step c hanges the velocit y distrib ution and, in this simplifi ed case, the space homogeneous problem admits the exa ct solution at time t + ∆ t f ( t + ∆ t ) = e − ∆ t/ε f ( t ) + (1 − e − ∆ t/ε ) E [ U ]( t ) . (38) The relaxation step of a Mont e Carlo metho d for t h e BGK equation consists in r ep lacing randomly selected particles with Maxwell ian particles with p robabilit y (1 − e − ∆ t/ε ). Thus v i ( t + ∆ t ) =  v i ( t ) , with pr ob ab ility e − ∆ t/ε E [ U ]( v ) , with probab ility 1 − e − ∆ t/ε (39) where E [ U ]( v ) represents a particle samp led from the Maxw ellian distribution with mo- men ts U . Th u s, fin ally , at eac h time step the m omen t guid ed Monte Carlo metho d reads as follo ws: (i) transp ort and co llide particles (37-39); (ii) solv e the fir st three momen t equ ations (24 ) and the additional equation for the third order moment (25); (iii) matc h the computed mass, m omen tum and energy of the particle solution (section 3.2) to those compu ted with th e momen t equ ations. Momen ts are reconstructed by sim p le summation form u las in eac h cell; fluxes are then obtained by in terp olat ion on the grid p oin ts and then discretized with Lax-F riedrichs t yp e cen tral sc hemes of fi rst or s econd ord er as describ ed in the previous sections. Remark 1 10 • Se c ond or der metho ds have b e en use d for the So d tests while first or der metho ds have b e en use d for al l others tests. We p oint out that se c ond or der scheme may pr o duc e lar ger fluctuations, esp e cial ly when slop e limiters ar e use d and that the switching te chnique b etwe en first and se c ond or der schemes as describ e d in (23) pr events the onset of these oscil lations in the c onsider e d test c ases. • After the r elaxation step (38), the p erturb ation term c an b e r ewritten as g ( t + ∆ t ) = f ( t + ∆ t ) − E [ U ( t + ∆ t )] = e − ∆ t/ε f ( t ) + (1 − e − ∆ t/ε ) E [ U ( t )] + − E [ U ( t + ∆ t )] = e − ∆ t/ε ( f ( t ) − E [ U ( t )]) = e − ∆ t/ε g ( t ) , (40) sinc e U ( t + ∆ t ) = U ( t ) i n the sp ac e homo gene ous c ase. Thus the discr etize d mo ment e qu ations (24) c an b e r ewritten as U n +1 i = U n i − ∆ t ∆ x ( ψ i +1 / 2 ( U n ) − ψ i − 1 / 2 ( U n )) + − ∆ t ∆ x e − ∆ t/ε Ψ i +1 / 2 ( < v m ( f n − E [ U n ]) > ) + + ∆ t ∆ x e − ∆ t/ε Ψ i − 1 / 2 ( < v m ( f n − E [ U n ]) > ) . (41) As ∆ t/ε gr ows, which me ans that the system appr o aches the e qu i librium, the c on- tribution of the k inetic term vanishes even though it is evaluate d thr ough p articles. This do es not happ en if we just c ompute the kinetic term ∂ x h v mg i fr om the p articles without c onsidering the structur e of the distribution fu nction f . This dr amatic al ly de cr e ases fluctuations when the Knudsen numb er is smal l. Sinc e this pr op erty is r elate d to the BGK structur e and we aim at a metho d that c an b e applie d to the ful l Boltzm ann e quation we do not take adva ntage of it in the numeric al r esults. We le ave the p ossibility to extend this ide a to the Boltzmann e qu ation using time r elaxe d Monte Carlo (TRMC) meth o ds[20] to futur e investigations. Accuracy t est First we rep ort on the resu lts of a sto chastic error analysis with resp ect to the num b er of particles. As reference solution w e co n sidered the a v erage of M in dep end en t realizations U M C = 1 M M X i =1 U i,M C (42) and U M G = 1 M M X i =1 U i,M G (43) where the tw o sub scripts M C an d M G indicate resp ectiv ely the reference solution for the Mon te Carlo metho d and for the Momen t Guided metho d. W e use t wo differen t reference 11 solutions since the tw o sc hemes presen t differen t discretization errors and th us they co n- v erge, w hen the num b er of particles go es to infinity , to different discretized solutions. Th e t w o r eference solutions are obtained by fixin g the time step and mesh size and le tting the n umb er of particles go to infinit y . In this wa y , b oth reference solutions con tain negligible sto c hastic error. A t the same time, b oth solutions in volv e space and time discretization errors. How ev er, the amoun t of su c h er r ors do es n ot change when the num b er of p articles v aries. Th erefore by comparing solutions obtained w ith a given ∆ t , ∆ x , but with a finite n umb er of particles to r eference solutions obtained with the same ∆ t , ∆ x , but with a v ery large num b er of particle, we obtain a true mea su re of the e r r or originating from the sto c hastic nature of the method . Th en, w e measure the quan tit y Σ 2 ( N ) = 1 M M X i =1 j max X j =1 ( U i,j − U j ) 2 (44) where U j represent s the reference s olution and j max represent s the num b er of mesh p oin t. The test c onsists of the follo win g initial data  ( x, 0) = 1 + a  sin 2 π x L u ( x, 0) = 1 . 5 + a u sin 2 π x L (45) 1 2 Z f | v | 2 dv = W ( x, 0) = 2 . 5 + a W sin 2 π x L with a  = 0 . 3 a u = 0 . 1 a W = 1 . This test p roblem giv es rise to a p erio d ic smo oth solution in the interv al t ∈ [0 , 5 × 10 − 2 ]. The results of this test in log-log scale are s h o wn in Figure 1. On the left, we rep orted the sto chastic error for the p ure Mon te Carlo and on the right for the Momen t Guided metho d. F r om top to b ott om, the errors for the three macroscopic quantit ies are depicted for different v alues of the Knudsen num b er. F or the Mon te Carlo case, the sto chastic error do es n ot s ubstanti ally c hange with resp ect to the K nudsen num b er and sho ws a con v ergence rate appr o ximativ ely equal to 1 / 2. At v ariance, for the Momen t Guid ed metho d, err ors decrease as the Knudsen num b er dimin ish and the con verge nce r ate of the metho d increases ac hieving v alues close to one. Th is b eha vior is due to th e fact that, for large Knudsen num b ers , the kin etic part of the solution, g , is not negligible and ev aluated through the DSMC method . By con trast, close to thermo d ynamical equilibrium, g → 0, whic h m eans that the Mon te Carlo comp onent of the solution carries only flu ctuations but no in f ormation. It is remark able that , in all analyzed regimes, the sto c hastic error of the Momen t Guided metho d is smaller than that of the pure p article solver. Unsteady sho c k test Next we consider an unsteady sho ck test ca se. Th is c hoice reflects the fact that the metho d is sp ecifically aimed at situations in wh ic h th e classical v ariance reduction tec hnique using 12 10 2 10 3 10 −4 10 −3 10 −2 10 −1 10 0 10 1 Density statistical error for MC method particles number N Σ ( ρ ) ε =1e−1 ε =1e−2 ε =1e−3 ε =1e −4 10 2 10 3 10 −4 10 −3 10 −2 10 −1 10 0 10 1 Density statistical error for HG method particles number N Σ ( ρ ) ε =1e−1 ε =1e−2 ε =1e−3 ε =1e −4 10 2 10 3 10 −3 10 −2 10 −1 10 0 10 1 Velocity statistical error for MC method particles number N Σ (u) ε =1e−1 ε =1e−2 ε =1e−3 ε =1e −4 10 2 10 3 10 −3 10 −2 10 −1 10 0 10 1 Velocity statistical error for HG method particles number N Σ (u) ε =1e−1 ε =1e−2 ε =1e−3 ε =1e −4 10 2 10 3 10 −2 10 −1 10 0 10 1 Temperature statistical error for MC method particles number N Σ (t) ε =1e−1 ε =1e−2 ε =1e−3 ε =1e −4 10 2 10 3 10 −2 10 −1 10 0 10 1 Temperature statistical error for HG method particles number N Σ (t) ε =1e−1 ε =1e−2 ε =1e−3 ε =1e −4 Figure 1: Statistica l error test: Solution at t = 0 . 05 for d ensit y (top), velocit y (middle) an d temp erature (b ottom). MC metho d (left), Hyd ro Guided MC metho d (r ight). Kn ud sen n umb er v ary from ε = 10 − 1 to ε = 10 − 4 . Squ ares indicate errors for ε = 10 − 1 , diamonds for ε = 10 − 2 , circles for ε = 10 − 3 while crosses indicate errors for ε = 10 − 4 . 13 time a ve r aging cannot b e used or tur ns out to b e u s eless, since time-a verag ing or using more p articles lea ds to the same co mp utational effort. Figures 2 to 5 consider the same initial data for the densit y , mean v elo cit y and temp er- ature with d ifferen t initial Knudsen n umber v alues, ranging from ε = 10 − 4 to ε = 10 − 1 . 100 particles p er cell are us ed and solutions are a v eraged o ver t wo different reali zations. Eac h Figure depicts the dens ity , mean v elo cit y and temp erature f rom top to b ottom, with the pu re Mon te Carlo solv er (on the left) and Moment Guided metho d (on the r igh t). In add ition, w e represent solutions of the compressib le Eu ler equations and as reference solution we used a discrete v elo cit y mo d el for the BGK equation [17]. These Figures sho w a large r eduction of fluctuations esp ecially for small Kn u dsen num b ers. So d sho c k tub e Finally w e lo ok at the classical So d Sh o c k T ub e test. F or this test case, we co n s ider the p ossibilit y of using second order flu id solv ers. W e observe that this c hoice has the effect of incr easing fluctuations far from thermo dynamical equilibrium . This is natural since w e miss the smo othing effect of a first order scheme. Ho wev er solutions obtained w ith firs t order sc hemes can b e un satisfactory in some situations b ecause of their large n umerical diffusion esp ecially close to thermod ynamical equilibrium. The solution whic h is adopted here consists in switc hing f r om the fir s t ord er to the second order s cheme according to the ratio of the th ermo dynamical flux w ith resp ect to the non equ ilibrium flu x. In pr actice, in eac h cell, the sc heme automatically u ses a sec ond ord er MUSCL scheme when th e kinetic term is small and a fi rst order sc h eme otherwise. Figures 6 to 9 consider the same initial data for the d ensit y , mean ve lo cit y an d tem- p erature with differen t initial K n ud sen n umb er v alues, whic h range from ε = 10 − 4 to ε = 10 − 1 . 100 p er cell are us ed and only one r ealizat ion is considered. As for the unsteady sho c k test eac h fi gure depicts the d ensit y , mean velocit y and temp erature from top to b ottom, with th e p u re Monte Carlo solv er (left) and the Momen t Guided m etho d (righ t). A reference solution obtained thr ough a d iscrete velocit y sc h eme [17] is represen ted in eac h figure as wel l as the solution of the compressible Eu ler equations. The figures sh o ws goo d results for a ll ranges of Knudsen n umb ers in terms of reduction of flu ctuations. Th e high order solve r do es n ot seem to increase th e v ariance bu t improv es the solution in the fluid limit. 5 Conclusions W e ha v e dev elop ed a new class of hybrid metho d s whic h aim at reducing the v ariance in Mon te Carlo sc hemes. The k ey idea consists in driving p article p ositions and velocities in suc h a wa y that momen ts giv en by the solution of the kinetic equatio n exactly matc h momen ts giv en by the solution of an app r opriate set of moment equations. It is imp ortan t to p oin t out th at the sc hemes wh ic h can be deriv ed through th is tec hn ique can b e easily implemen ted in existing Mon te Carlo co d es through few mo difi cations: adding a fluid solv er a n d a routine for the momen t matc hin g. 14 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.0001 x ρ (x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.0001 x ρ (x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.0001 x u(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.0001 x u(x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.0001 x T(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.0001 x T(x,t) MG DVM rif EULER rif Figure 2: Unsteady Sho ck: Solution at t = 0 . 065 for the density (top), v elo cit y (mid- dle) and temp erature (b ottom). MC metho d (le ft), Momen t Guided method MG (righ t). Kn u dsen num b er ε = 10 − 4 . Reference solution: d ash dotte d line. Eu ler solutio n : con tin- uous line. Monte Carlo or Mo ment Guided: circles p lus con tinuous line. 15 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.001 x ρ (x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.001 x ρ (x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.001 x u(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.001 x u(x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.001 x T(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.001 x T(x,t) MG DVM rif EULER rif Figure 3: Unsteady Sho ck: Solution at t = 0 . 065 for the density (top), v elo cit y (mid- dle) and temp erature (b ottom). MC metho d (le ft), Momen t Guided method MG (righ t). Kn u dsen num b er ε = 10 − 3 . Reference solution: d ash dotte d line. Eu ler solutio n : con tin- uous line. Monte Carlo or Mo ment Guided: circles p lus con tinuous line. 16 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.01 x ρ (x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.01 x ρ (x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.01 x u(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.01 x u(x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.01 x T(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.01 x T(x,t) MG DVM rif EULER rif Figure 4: Unsteady Sho ck: Solution at t = 0 . 065 for the density (top), v elo cit y (mid- dle) and temp erature (b ottom). MC metho d (le ft), Momen t Guided method MG (righ t). Kn u dsen num b er ε = 10 − 2 . Reference solution: d ash dotte d line. Eu ler solutio n : con tin- uous line. Monte Carlo or Mo ment Guided: circles p lus con tinuous line. 17 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.1 x ρ (x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 Density, Knudsen number=0.1 x ρ (x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.1 x u(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 Velocity, Knudsen number=0.1 x u(x,t) MG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.1 x T(x,t) MC DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 3.5 4 4.5 5 5.5 6 6.5 7 7.5 Temperature, Knudsen number=0.1 x T(x,t) MG DVM rif EULER rif Figure 5: Unsteady Sho ck: Solution at t = 0 . 065 for the density (top), v elo cit y (mid- dle) and temp erature (b ottom). MC metho d (le ft), Momen t Guided method MG (righ t). Kn u dsen num b er ε = 10 − 1 . Reference solution: d ash dotte d line. Eu ler solutio n : con tin- uous line. Monte Carlo or Mo ment Guided: circles p lus con tinuous line. 18 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.0001 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.0001 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.0001 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.0001 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.0001 x T(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.0001 x T(x,t) HG DVM rif EULER rif Figure 6: Sod S ho c k T ub e T est: Solutio n at t = 0 . 05 f or the dens ity (top), ve lo city (middle) and temp erature (b ottom). MC m etho d (left), Momen t Guided MG metho d (righ t). Knudsen n u m b er ε = 10 − 4 . Reference solution: dash dotted line. Euler sol u tion: con tin uous line. Mon te Carlo or Momen t Guided: circles plus con tinuous line. 19 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.001 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.001 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.001 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.001 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.001 x T(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.001 x T(x,t) HG DVM rif EULER rif Figure 7: Sod S ho c k T ub e T est: Solutio n at t = 0 . 05 f or the dens ity (top), ve lo city (middle) and temp erature (b ottom). MC m etho d (left), Momen t Guided MG metho d (righ t). Knudsen n u m b er ε = 10 − 3 . Reference solution: dash dotted line. Euler sol u tion: con tin uous line. Mon te Carlo or Momen t Guided: circles plus con tinuous line. 20 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.01 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.01 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.01 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.01 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.01 x T(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.01 x T(x,t) HG DVM rif EULER rif Figure 8: Sod S ho c k T ub e T est: Solutio n at t = 0 . 05 f or the dens ity (top), ve lo city (middle) and temp erature (b ottom). MC m etho d (left), Momen t Guided MG metho d (righ t). Knudsen n u m b er ε = 10 − 2 . Reference solution: dash dotted line. Euler sol u tion: con tin uous line. Mon te Carlo or Momen t Guided: circles plus con tinuous line. 21 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.1 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 1.2 Density, Knudsen number=0.1 x ρ (x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.1 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1 −0.5 0 0.5 1 1.5 2 2.5 Velocity, Knudsen number=0.1 x u(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.1 x T(x,t) HG DVM rif EULER rif 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 2 3 4 5 6 7 8 9 10 11 Temperature, Knudsen number=0.1 x T(x,t) HG DVM rif EULER rif Figure 9: Sod S ho c k T ub e T est: Solutio n at t = 0 . 05 f or the dens ity (top), ve lo city (middle) and temp erature (b ottom). MC m etho d (left), Momen t Guided MG metho d (righ t). Knudsen n u m b er ε = 10 − 1 . Reference solution: dash dotted line. Euler sol u tion: con tin uous line. Mon te Carlo or Momen t Guided: circles plus con tinuous line. 22 Preliminary n u merical results sho w redu ctions of flu ctuations in all regimes compared to DSMC. The redu ction b ecomes stronger as we approac h equilibriu m. Numerical con- v ergence tests sho w b etter p erformances of the prop osed method , in terms of sto c hastic error, compared to pure Monte Carlo sc hemes. F or these problems the momen t guided metho d seems very promising, leading to solutions w h ic h conta in less fl u ctuations at a computational cost whic h is comparable to the cost of a traditional Mon te Carlo solv er, and the ad d ition of the cost of a macroscopic solv er f or the compressible Euler equ ations, whic h is usually computationally less exp ensive than the Monte Carlo metho d . Currently , we are working on extensions of the pr esen t method to the full Boltzmann equation in the multidimensional case. T o this aim we plan to u se b oth cla ssical Monte Carlo metho ds lik e Bird or Nan bu metho ds [18] and time relaxed Mon te Carlo (TRMC) tec hniques [20]. Moreov er w e plan to explore other p ossible algorithms w h ic h can p ossibly further reduce fluctuations, suc h as matc hing higher order momen ts and/or using higher order closure of the h ierarc hy in order to solv e a larger s et of h ydr o dynamics equations, or using hybrid r epresen tations of the d istribution function [5]. W e hop e to b e able to present other results su pp orting this metho d ology in the near future [7]. References [1] G.A.Bird , Mole cular gas dynamics and dir e ct simulation of gas flows , Clare ndo n P ress, Oxford (199 4). [2] C.K. Birsd all, A .B. Langdon , Plasma Physics Via Computer Simulation , Institute of Physics (IOP), Series in Plasma P hysics (2004 ). [3] J. Bur t, I . Boyd , A low diffusion p article metho d for simulating c ompr essible inviscid flows , Journal of Co mputational Physics, V olume 2 27, Issue 9, 20 April 2008 , pp. 4653- 4670 [4] R. E. Caflisch , Monte Carlo and Quasi-Monte Carlo Metho ds , Acta Numerica (19 98) pp. 1 – 49. [5] R. E. Caflisch, L. P areschi , T owar ds an hybrid metho d for r ar efie d gas dynamics , IMA V ol. App. Math., vol. 135 (200 4), pp. 57–73 . [6] C. Cercignani , The Boltzmann Equation and Its Applic ations , Springe r -V erlag, New Y o rk, (1988). [7] P. Deg ond , G. Dimarco, L. P areschi , Moment Guide d Monte Carlo Schemes for the Boltzmann e quation , W ork in Progr ess 2009 . [8] P.Degond, J. -G. Liu, L. Mieussens , Macr osc opic fluid mo dels with lo c alize d kinetic up- sc aling effe cts , SIAM MMS, vol. 5 (2006 ), pp. 940- 979. [9] G. Dimarco, L . P areschi , A Fluid Solver In dep endent Hybrid metho d for Mult isc ale Kinetic Equations , SIAM J. Sci. Comput. (Submitted). [10] H.Grad , On the kinetic the ory of ra r efie d gases , Co mmun. Pure Appl. Ma th. 2 , 33 1 (19 49). [11] W. E, B. En gquist , The heter o gene ous multisc ale met ho ds , Comm. Math. Sci., vol. 1 (20 03), pp. 87-1 33. 23 [12] T. Hom olle, N. Hadjiconst antinou , A low-varianc e deviational simulation Monte Carlo for the Boltzmann e quation. J ournal of Computatio nal Physics 226 (200 7), pp 2341 –2358 . [13] T. Homoll e, N. Hadjiconst antinou , L ow-varianc e devi ational simulation Monte Carlo. Physics of Fluids 19 (2007) 04 1701 . [14] D. Levermore , Moment clo su r e hier ar chies for kinetic t he ories , Jo urnal of Statistica l Physics, V o l. 83 , No. 5-6, (199 6). [15] R. J. L eVeque , Numeric al Metho ds for Conservation L aws , Lectures in Mathematics, Birkhauser V erlag , Basel (1 9 92). [16] S. Liu , Monte Ca rlo strateg ies in scientific computing, Springer, (2004 ). [17] L. Mieussens , Discr ete V elo city Mo del and Impli cit Scheme for t he BGK Equat ion of Ra r efie d Gas Dynamic ,Mathematical Mo dels and Me tho ds in Applied Scienc e s , V ol. 10 No. 8 (2 0 00), 1121– 1149 . [18] K. Nanbu , Dir e ct simulation scheme derive d fr om the Boltzmann e quation , J. Phys. So c. Japan, vol. 49 (1 9 80), pp. 204 2 –204 9. [19] I. M ¨ uller, T. Rug geri , Ra t ional Extende d Thermo dynamics. Springer T racts in Natural Philosophy V ol. 3 7 pp. 3 93, I I Edition - Springer V erla g (1 9 98). [20] L. P areschi, G. Russo , Time R elaxe d Monte Carlo metho ds for the Boltzmann e quation , SIAM J. Sci. Co mput. 23 (200 1), pp. 1253 –1273 . [21] L. P areschi , Hybrid multisc ale metho ds for hyp erb olic and kinetic pr oblems , Esaim Pro ceed- ings, V o l. 15 , T. Goudon, E. So nnendruck er & D. T alay Editor s, pp.87-12 0, (200 5 ). [22] D. I. Pullin , Dir e ct simulation metho ds for c ompr essible inviscid ide al gas flow , J. Co mput. Phys., 34 (198 0), pp. 231 –244 . 24