Overview Of Accelerated Simulation Methods For Plasma Kinetics

Overview of Accelerated Simulation Methods forPlasma KineticsR.E. Caflisch1In collaboration with: J.L. Cambier2 , B.I. Cohen3 , A.M.Dimits3 , L.F. Ricketson1,4 , M.S. Rosin1,5 , B. Yann11 UCLAMath DeptInstitute4 Courant2 AFRL5 Pratt3 LLNLInstituteAFRL Program ReviewJanuary 20, 2015Work performed under funding from AFOSR and DOE.Russel CaflischAccelerated Simulation for Plasma Kinetics

IntroductionRussel CaflischAccelerated Simulation for Plasma Kinetics

Kinetic TheoryPhase space particle number density f ( x , v , t)Boltzmann Equation: t f v · x f a · v f C (f , f )(1)C (f , f ) is:Boltzmann collision operator for rarefied gas dynamics (RGD)Landau-Fokker-Planck operator for Coulomb collisionsOperator for excitation/deexcitation, ionization/recombinationfor collisional/radiative (CR) kineticsThis talk omits spatial variation and mean fields.Russel CaflischAccelerated Simulation for Plasma Kinetics

Binary Collision MethodCollisions are simulated directly as in DSMC, withDistribution function is represented by set of discrete particlesFor short range interactions (RGD and CR kinetics), eachparticle has ν t collisions in time step t for collision rate ν.Collision partners and parameters are randomly chosen.For long range interactions (Coulomb), each numericalcollision is an aggregate of grazing collisions over time step.Every particle collides once in every time stepCollisions are aggregates depending on tRussel CaflischAccelerated Simulation for Plasma Kinetics

Computational Obstacles for Monte Carlo SimulationComputational complexity for high collisionality (nearcontinuum limit)Hybrid method for RGD and Coulomb collisionsNegative particle methodMany time stepsMulti-Level Monte Carlo (MLMC)Multiscale featuresCR kineticsRussel CaflischAccelerated Simulation for Plasma Kinetics

Hybrid Schemes for Binary CollisionMethodsRussel CaflischAccelerated Simulation for Plasma Kinetics

Hybrid SchemeCombine fluid and particle simulation methods1 :Treat as particlesSeparate f into Maxwellianand non-Maxwelliancomponents: f m fpTreat m as fluid solvesEuler equationsSimulate fp by Monte CarloalgorithmTreat as fluidInteraction of m and fp :sample particles from m andcollide them with particlesfrom fp .1Caflisch et. al, Multiscale Model. Simul. 7 (2008) 865-887.Russel CaflischAccelerated Simulation for Plasma Kinetics

Thermalization/dethermalizationHybrid scheme is efficient because collisions between m and mneed not be simulated. Efficiency increased by thermalizationCollisions drive particles into equilibriumMove particles from fp to m when they have collided enough(thermalization)Move sampled particles from m into fp if the collision isstrong enough (dethermalization)(De)Thermalization criterion using entropy2Alternative criterion based on scattering angle3Generalization to spatial inhomogeneities is nontrivial23Ricketson et. al, J. Comp. Phys., 2014.Dimits et. al., private communicationRussel CaflischAccelerated Simulation for Plasma Kinetics

Bump-on-Tailt 1.2tFP0.20.150.150.150 10.1f0.10.050.05 0.50vx0.50 110.10.05 0.5t 3.6tFP0vx0.50 110.20.150.150.150.1f0.10.050.05 0.50vx0.510 10vx0.510.51t 11tFP0.20 1 0.5t 6tFP0.2fft 2.4tFP0.2fft 00.20.10.05 0.5Russel Caflisch0vx0.510 1 0.50vxAccelerated Simulation for Plasma Kinetics

Scheme ComparisonBump on Tail Test Problem0.060.055Relative L1 Error0.050.0450.040.0350.03Entropy SchemeEntropy DethermalizationScattering Angle SchemeScattering Angle Dethermalization0.0250.025101520253035Improvement Factor vs. PICRussel CaflischAccelerated Simulation for Plasma Kinetics40

Negative ParticlesRussel CaflischAccelerated Simulation for Plasma Kinetics

Hybrid methods need negative particlesThe hybrid method presented above, assumes that f m, This may beinefficient if there is a “defect”’ in the Maxwellian.1f0.80.60.40.20 6 4 2024681012vRussel CaflischAccelerated Simulation for Plasma Kinetics

Hybrid methods need negative particlesThe hybrid method presented above, assumes that f m, This may beinefficient if there is a “defect”’ in the Maxwellian.If f m fp (with fp 0 so that fp can represented by particles) then m willbe small.Russel CaflischAccelerated Simulation for Plasma Kinetics

Hybrid methods need negative particlesThe hybrid method presented above, assumes that f m, This may beinefficient if there is a “defect”’ in the Maxwellian.Better:f m fp fnwith fp 0, fn 0. Represent the defect by “negative particles”!Russel CaflischAccelerated Simulation for Plasma Kinetics

Meaning of Negative ParticlesA negative particle w (in fm ) cancels a particle in m (or infp ) which should not be present.A collision P N between a positive particle v and anegative particle w cancels a collision P P between v anda positive particle w that should not have occurred.The P P collision removes v and w and adds v 0 and w 0 :P-P:v , w v 0 , w 0So the P N collision should add a v , remove w (equivalentto adding w ) and add negative particles to cancel v 0 and w 0P-N:00v , w 2v , v , w This can be derived from the Boltzmann equation(Hadjiconstantinou 2005).Russel CaflischAccelerated Simulation for Plasma Kinetics

General rules for collisions with negative particlesP-P: v , w v 0 , w 0P-N:N-N:P-M:v , w 2v , v 0 , w 0v , w 2v , 2w , v 0 , w 0m, v m, w , v 0 , w 0N-M: m, v m, w , v 0 , w 0 .Particle number increases!Since every particle collides in every time step, it’s muchworse for Coulomb collisions!Bokai Yan will present a new method for overcoming thisparticle increase.Russel CaflischAccelerated Simulation for Plasma Kinetics

Multilevel Monte Carlo (MLMC) forLangevin MethodRussel CaflischAccelerated Simulation for Plasma Kinetics

Langevin FormulationLinear LFP equation for f (v , t) is equivalent to stochasticdifferential equations (SDEs) for v (t)dvi Fi dt Dij dWj ,(2)where f is probability density of v and i, j are component indicesW W (t) is Brownian motion in velocitydW is white noise in velocityDirect extension to spatial dependenceValid for nonlinear LFP, if F and D are updated as neededRussel CaflischAccelerated Simulation for Plasma Kinetics

Discretization of SDEsEuler-Maruyama discretization in time:vi,n 1 vi,n Fi,n t Dij,n Wj,n ,(3) Wn Wn 1 Wn(4)in which vi,n vi (tn ) and Fn F(vn )Computational cost vs. error ε:Statistical error is O(N 1/2 ) t error is O( t), since W O( t) and randomOptimal choice is ε N 1/2 tCost N t 1 ε 3Russel CaflischAccelerated Simulation for Plasma Kinetics

MLMC scheme combines multiple solutions with varyingt, NlMLMCBasicsParametric IntegrationConvergent sum, L ! 1Like multi-gridLXwhen using Milstein methodNLv̂L E[v0 ] E[(vl vl 1 )]This can be repeated onl 1multiple levels (perhaps usingIntroducestep levels, t T 2 , for 0,smooth)., LhigherordertimeinterpolationExpectationswith Nl samples if f (x, λ) is sufficientlyStatistical error converges with tmethod.!Grid!!"!!!!!!!!!!!!!!!!!! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !"λ7Giles in “Monte Carlo and Quasi-Monte Carlo Method”, Springer-Verlag, (2006)Use1. computationat coarse level as a “control variate” forThis doesn’t quite fit into the multilevel framework I’vecomputationleveldescribed,but atthefinercomplexityanalysis is very similar.Optimally combine the different levelsMultilevel Monte Carlo – p. 10Motivated by mutigrid methodsRussel CaflischAccelerated Simulation for Plasma Kinetics

MLMC ScalingFor RMSE ε, the complexity now scales like4 O ε 2 (log ε)2 for Euler-MaruyamaCost O ε 2for Milstein(5)Notes:MLMC-Euler-Maruyama scales better than standard MCMLMC-Milstein is even better O ε 2 scaling is possible without Milstein, using antitheticsampling method 545Giles, Operations Research, 56(3):607, 2008Giles & Szpruch, arXiv:1202.6283, 2012Russel CaflischAccelerated Simulation for Plasma Kinetics

A Sample Plasma ProblemDirectMLMC EulerMLMC Milstein (lnε)2K ε210010-110-510-410-3εRosin, Ricketson, et. al., J Comp Phys 2014Russel CaflischAccelerated Simulation for Plasma Kinetics

Direct Simulation of CR KineticsRussel CaflischAccelerated Simulation for Plasma Kinetics

CR KineticsKinetic simulation of electron-impact excitation/deexcitation andionization/recombination.6Plasma representationDiscrete particles for electron energy distribution f (E )Continuum densities ρk for atoms at k m 1 levels ofexcitation and ρi for ions (Bohr model)Cross-sectionsSemi-classical cross-sections for excitation and ionization.Detailed balance yields cross-sections for deexcitation andrecombination.Closely related to work of Hai Le.6Yan et. al., J Comp Phys to appearRussel CaflischAccelerated Simulation for Plasma Kinetics

Boltzmann EquationElectronsee t f (E ) QIR QEDAtomic level kkk t ρk QIR QEDIonsi t ρi QIRQ is the collision operator for ionization/recombination (IR) andexcitation/deexcitation (ED).Russel CaflischAccelerated Simulation for Plasma Kinetics

Entropy InequalityEntropy S isS kN(log (N/G ) 1)N number density of particlesG degeneracy of particlesk Boltzmann constantBoltzmann H-theorem, for H -S, isdH 0dtwith equality only for equilibriumNew(?) explicit formula for S in terms of distributionfunctions.Multiple non-Maxwellian equilibria for excitation/deexcitation!H theorem is always valid.Russel CaflischAccelerated Simulation for Plasma Kinetics

Multiscale and Singular FeaturesMultiscale and singular features are a bottleneck in Monte CarlosimulationMany processesm states (m 1 levels ion)interactions between any two states k, k 0 (including k i)m(m 1)/2 different types of interactionsWide range of interaction ratesinteraction k (k 1) has rate O(k 6 )!Singularity in recombination rate for recombining electronenergies E1 and E2 1rrec O E1 E2Russel CaflischAccelerated Simulation for Plasma Kinetics

Numerical ResultsEvolution for ionization/recombination only27t 0.00 ps845atom distributionx 104actualequilibrium6410123Levels45027t 0.02 ps00.51Energy / IH1.520.51Energy / IH1.520.51Energy / IH1.5244x 10153x 10102510123Levels450274t 0.50 psactualequilibrium224electron distributionx 10323x 101.52110.50045x 10123Levels4500Maxwellian equilibrium as expectedSimilar results if excitation/deexcitation is includedRussel CaflischAccelerated Simulation for Plasma Kinetics

Numerical ResultsEvolution for ionization/recombination onlyentropy010 110 210 310 410 51000.050.10.150.20.25time / ps0.30.350.40.450.5Decrease in H as expectedRussel CaflischAccelerated Simulation for Plasma Kinetics

Numerical ResultsEvolution for excitation/deexcitation only, with 2 atomic levels26t 0.00 ps846atom distributionx 10actualequilibrium261.5412electron distributionx 10actualequilibrium0.5012000.511.52Energy / IH2.533.50.511.52Energy / IH2.533.50.511.52Energy / IH2.533.5Levels26t 0.02 ps845x 10156x 101045201200Levels26t 0.52 ps645x 101541020x 1051200LevelsEquilibrium is clearly non-MaxwellianFor 10 levels, equilibrium very close to MaxwellianRussel CaflischAccelerated Simulation for Plasma Kinetics

Numerical ResultsEvolution for excitation/deexcitation only, with 2 atomic levelsEntropy010 110 210 31000.10.20.30.40.50.60.7time / psDecrease in H as expected, in spite of non-MaxwellianRussel CaflischAccelerated Simulation for Plasma Kinetics

Summary: Accelerated Methods for Monte CarloSimulationHybrid method for RGD and Coulomb collisionsNegative particle methodMulti-Level Monte Carlo (MLMC)CR kineticsRussel CaflischAccelerated Simulation for Plasma Kinetics

Boltzmann collision operator for rare ed gas dynamics (RGD) Landau-Fokker-Planck operator for Coulomb collisions Operator for excitation/deexcitation, ionization/recombination for collisional/radiative (CR) kinetics This talk omits spatial variation and mean elds. Russel Ca isch Accelerated Simulation for Plasma Kinetics