TN-02_Numerical1DDriftKineticPeriodicBoundaryConditions
=======================================================

.. meta::
   :description: technical note

   :keywords: Report,2047357-TN-02,,M2.1,Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,M.,Barnes1,,F.,I.,Parra1,and,M.,R.,Hardman1,1,Rudolf,Peierls,Centre,for,Theoretical,Physics,,University,of,Oxford,,Clarendon,Laboratory,,Parks,Road,,Oxford,OX1,3PU,,UK,E-mail:,michael.barnes@physics.ox.ac.uk,1.,Introduction,We,expect,that,one,of,the,biggest,challenges,in,numerically,solving,drift,kinetic,equations,in,the,plasma,edge,is,treating,the,motion,of,electrons,along,the,magnetic,ﬁeld.,Because,the,electrons,are,light,,they,move,rapidly,along,the,ﬁeld,,placing,a,severe,stability,restriction,on,the,step,size,for,explicit,time,advance,schemes.,Unfortunately,,an,implicit,treatment,is,not,straightforward,due,to,an,implicit,dependence,of,the,electrostatic,potential,on,the,charged,particle,distribution,functions.,One,of,the,main,aims,of,our,research,is,to,develop,and,test,a,novel,analytical,model,and,associated,numerical,algorithm,for,relaxing,this,restriction.,As,a,ﬁrst,step,towards,this,goal,,we,have,developed,a,new,code,in,the,programming,language,Julia,to,simulate,a,simple,model,for,parallel,dynamics,(described,in,our,Jan,2021,report,[1]),without,the,novel,moment-based,approach,that,we,intend,to,ultimately,employ.,This,code,will,be,used,to,test,diﬀerent,numerical,approaches,and,will,be,built,upon,to,create,the,code,for,the,moment-based,approach.,It,will,also,be,a,useful,benchmark,against,which,the,moment-based,approach,can,be,tested.,In,this,report,we,give,a,brief,description,of,a,simple,model,for,parallel,plasma,dynamics,before,describing,the,corresponding,code,and,presenting,a,comparison,of,numerical,results,with,the,analytical,benchmark,we,have,derived.,2.,Model,equations,A,detailed,derivation,of,the,model,we,consider,is,provided,in,the,previous,(Jan,2021),report,[1].,Here,we,provide,a,brief,overview,of,the,model,for,the,Reader’s,convenience.,The,model,we,consider,consists,of,a,single,ion,species,of,charge,e,,a,single,neutral,species,,and,an,electron,species,modelled,as,having,a,Boltzmann,response,,all,immersed,in,a,straight,,uniform,magnetic,ﬁeld,in,the,z,direction.,We,allow,for,charge,exchange,collisions,between,ions,and,neutrals,but,do,not,account,for,intra-species,collisions.,Finally,,we,assume,that,the,plasma,is,homogeneous,in,the,plane,perpendicular,to,the,1,Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,2,magnetic,ﬁeld.,With,these,assumptions,,our,model,system,of,equations,is,and,.,=,∂fi,∂t,+,v,(cid:107),∂fi,∂z,−,e,mi,∂φ,∂z,∂fi,∂v,(cid:107),=,Rin,(nnfi,−,−,nifn),,,∂fn,∂t,+,v,∂fn,∂z,(cid:107),=,ns(z,,t),=,−,∞,Rin,(nifn,−,nnfi),,,dv,(cid:107),fs(z,,v,,,t),,(cid:107),(cid:90),−∞,ni,=,Ne,exp,,,eφ,Te,(cid:19),(cid:18),(1),(2),(3),(4),⊥,v,⊥,dϑdv,and,v,with,fs,Fs,the,marginalized,particle,distribution,function,for,species,s,,v,the,components,of,the,particle,velocity,parallel,and,perpendicular,to,the,(cid:82),(cid:107),magnetic,ﬁeld,,respectively,,ϑ,the,gyro-angle,,mi,the,ion,mass,,t,the,time,,φ,the,electrostatic,potential,,and,Rin,the,charge,exchange,collision,frequency.,⊥,For,our,boundary,conditions,,we,impose,periodicity,on,fs,in,both,z,and,v,,,with,(cid:107),periods,Lz,and,Lv(cid:107),,respectively.,There,is,also,the,option,to,impose,zero,boundary,at,the,upwind,boundary,of,the,domain.,As,fs,should,go,to,conditions,on,z,and,v,(cid:107),zero,in,v,,,imposition,of,zero,boundary,conditions,and,periodic,boundary,as,v,(cid:107),conditions,should,be,equivalent,as,long,as,Lv(cid:107),is,suﬃciently,large.,Note,that,with,either,choice,of,boundary,conditions,,the,line-averaged,density,(cid:107),→,±∞,Lz,0,dz,ns,is,conserved.,(cid:82),,,We,normalize,Eqs.,(1)-(4),by,deﬁning,˜fs,.,=,fs,,,,,vth,i√π,Ne,vth,i,Lz,z,Lz,v,(cid:107),vth,i,ns,Ne,eφ,Te,,,,,,,.,=,t,˜t,.,=,˜z,.,=,˜v,(cid:107),.,=,.,=,˜ns,˜φ,and,with,vth,i,.,=,˜Rin,.,=,Rin,Ne,Lz,vth,i,2Te/mi.,In,terms,of,these,normalised,quantities,,Eqs,(1)-(4),become,(cid:112),∂,˜fi,∂˜t,∂,˜fi,∂,˜z,−,∂,˜φ,∂,˜z,1,2,+,˜v,(cid:107),∂,˜fi,∂˜v,(cid:107),=,˜Rin,−,˜nn,˜fi,−,(cid:16),˜ni,˜fn,,,(cid:17),(5),(6),(7),(8),(9),(10),(11),(12),Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,3,and,∂,˜fn,∂˜t,∂,˜fn,∂,˜z,+,˜v,(cid:107),=,e,˜φ,=,˜ni,=,˜Rin,−,1,√π,˜ni,(cid:16),∞,˜nn,˜fn,−,˜fi,,d˜v,(cid:107),(cid:90),−∞,˜nn,=,1,√π,∞,(cid:90),−∞,˜fn.,d˜v,(cid:107),˜fi,,,(cid:17),(13),(14),(15),3.,Numerical,implementation,The,algorithms,described,in,this,Section,have,been,implemented,in,the,code,,written,in,the,Julia,programming,language,,currently,available,on,GitHub,at,https://github.,com/mabarnes/moment_kinetics.,3.1.,Time,advance,We,evolve,Eqs.,(12)-(15),using,a,time-marching,scheme,(as,opposed,to,an,eigensolver),due,to,its,eﬃciency,and,to,the,nonlinear,nature,of,the,system,of,partial,diﬀerential,equations.,In,particular,,we,employ,a,member,of,the,family,of,Strong,Stability,Preserving,(SSP),Runge-Kutta,(RK),schemes;,see,,e.g.,,[2,,3,,4].,Current,SSPRK,options,implemented,in,the,code,are,SSPRK1,(forward,Euler),,SSPRK2,(Heun’s,method),and,SSPRK3,(Shu-Osher,method).,The,user,can,also,specify,the,use,of,‘ﬂip-,ﬂop’,Lie,operator,splitting,,described,in,Appendix,A.,Operator,splitting,limits,the,time,advance,scheme,to,second,order,accuracy,in,step,size,,but,could,be,useful,for,separately,treating,diﬀerent,pieces,of,physics.,Here,we,describe,the,current,default,option,,which,is,the,Shu-Osher,method,(SSPRK3),without,operator,splitting.,For,convenience,of,notation,,we,express,the,normalised,drift,kinetic,equations,for,the,ions,and,neutrals,in,the,vector,form,∂f,∂t,=,G[f,],,(16),with,f,=,(,˜fi,,˜fn)T,and,G,the,drift,kinetic,operator,account,for,parallel,streaming,,parallel,acceleration,(for,the,ions),and,charge,exchange,collisions.,The,Shu-Osher,method,for,advancing,this,system,of,equations,is,a,3-stage,,SSPRK,scheme,that,is,3rd,order,accurate,in,time,step,size,∆t.,It,is,given,by,f,(1),=,f,n,+,∆tG,[f,n],,,f,(2),=,f,n+1,=,3,4,1,3,f,n,+,f,n,+,1,4,2,3,f,(1),+,∆tG,f,(1),(cid:0),f,(2),+,∆tG,(cid:2),f,(2),(cid:3)(cid:1),,,,,(17),where,the,superscript,n,denotes,the,time,level.,(cid:0),(cid:2),(cid:3)(cid:1),Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,4,3.2.,Spatial,discretisation,There,are,two,discretisation,schemes,implemented,in,the,code:,ﬁnite,diﬀerences,and,Chebyshev,(pseudo)spectral,elements.,The,user,can,choose,at,run-time,which,scheme,to,use,for,each,of,the,z,and,v,(cid:107),coordinates.,3.2.1.,Finite,diﬀerence,discretisation.,For,the,ﬁnite,diﬀerence,discretisation,,the,corresponding,coordinate,grid,is,uniform,on,the,domain,[,L/2,,L/2],,with,L,the,coordinate,box,length.,The,default,method,employed,for,derivatives,is,3rd,order,upwind,diﬀerences,,though,1st,and,2nd,order,schemes,are,also,available,as,options.,For,an,overview,of,upwind,diﬀerences,and,a,discussion,of,the,merits,of,the,diﬀerent,upwind,schemes,,see,,e.g.,[5].,The,associated,integration,weights,used,for,ﬁeld-line,averages,in,z,integration,required,for,obtaining,ﬁelds/moments,are,obtained,using,and/or,for,the,v,(cid:107),the,composite,Simpson’s,rule,(sometimes,referred,to,as,composite,Simpson’s,1/3,rule:,−,L,0,(cid:90),dx,f,(x),h,3,≈,(N,−,1)/2,j=1,(cid:88),(f,(x2j,1),+,4f,(x2j),+,f,(x2j+1)),,,−,(18),1),is,the,where,N,is,the,number,of,grid,points,in,the,coordinate,x,,and,h,=,L/(N,uniform,grid,spacing.,The,composite,rule,(18),is,only,applicable,for,N,odd,,so,it,is,supplemented,at,the,boundary,by,Simpson’s,3/8,rule,when,N,is,even.,−,3.2.2.,Chebyshev,spectral,elements.,When,using,Chebyshev,spectral,elements,(see,,e.g.,,[6]),,the,corresponding,coordinate,grid,is,the,Gauss-Chebyshev-Lobatto,grid,on,each,element.,For,a,description,of,Chebyshev-Gauss,quadrature,,see,,e.g.,[7].,Inclusion,of,the,endpoints,within,each,element,facilitates,enforcement,of,continuity,at,element,boundaries,,and,the,use,of,Chebyshev,polynomials,as,a,basis,enables,the,use,of,Fast,Fourier,Transforms.,In,our,code,,these,transforms,are,done,using,the,widely-used,FFTW,library,[8].,The,associated,integration,weights,used,for,ﬁeld-line,averages,in,z,and/or,for,the,v,integration,required,for,obtaining,ﬁelds/moments,are,obtained,using,Clenshaw-,(cid:107),Curtis,quadrature,rules,[9].,Clenshaw-Curtis,quadrature,is,convenient,,as,it,allows,for,the,use,of,endpoints,in,the,integration,domain,(which,is,dictated,by,the,use,of,a,Gauss-,Chebyshev-Lobatto,grid),while,still,exactly,integrating,polynomials,up,to,degree,N,1,,with,N,the,number,of,points,within,the,element.,−,4.,Numerical,results,implementation,of,Eqs.,(12)-(15),,we,compare,our,To,benchmark,our,numerical,simulation,results,with,the,analytical,benchmarks,developed,in,[1].,In,particular,,we,have,initialised,the,distribution,functions,for,the,ions,and,neutrals,to,be,of,the,form,˜fs,=,ns,Ne,(cid:18),Te,T,s,(cid:19),1/2,exp,˜v2,(cid:107),−,(cid:18),Te,T,s,(cid:19),,,(19),Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,5,with,ns,=,ns,+δns,,and,an,overline,denoting,a,ﬁeld,line,average.,The,piece,of,the,density,that,varies,along,z,,δns,,is,chosen,to,be,small,compared,to,ns,(δns/ns,=,0.001),so,that,the,system,of,equations,can,be,linearised,to,a,good,approximation.,This,facilitates,comparisons,with,the,linear,analytical,theory,presented,in,[1].,For,all,cases,shown,here,,ni,=,nn,=,Ne/2,,T,i,=,T,n,and,mi,=,mn.,Both,the,electron-ion,temperature,ratio,and,the,charge,exchange,collision,frequency,are,varied,,and,damping,rates,and,frequencies,are,extracted,by,considering,the,time,evolution,of,the,spatially-varying,component,of,the,electrostatic,potential,,δφ.,In,particular,,a,least-squares,ﬁt,for,δφ(t)/δφ(t0),is,done,ϕ),for,each,simulation,to,a,function,of,the,form,exp(γ(t,to,obtain,the,damping,rate,γ,,frequency,ω,and,phase,ϕ.,The,results,are,given,in,Fig.,1.,There,is,excellent,agreement,across,a,wide,range,of,temperature,ratios,and,charge,exchange,collision,frequencies,,both,for,the,damping,of,ﬁnite,frequency,modes,(corresponding,to,the,solid,lines),and,to,a,zero,frequency,mode,that,appears,at,larger,collisionalities,(dashed-dotted,lines).,ϕ)/,cos(ωt0,t0)),cos(ωt,−,−,−,−,Figure,1.,Normalized,growth,rate,and,real,frequency,as,a,function,of,the,ion-electron,temperature,ratio.,Note,that,the,normalising,vth,employed,here,diﬀers,from,the,vth,i,employed,in,the,text:,It,is,chosen,to,be,vth,=,2Ti/mi,to,facilitate,comparison,with,the,analytical,results,obtained,in,[1].,(cid:112),The,minor,discrepancies,between,the,analytical,and,numerical,damping,rates,that,are,apparent,for,a,handful,of,the,cases,are,due,to,the,simultaneous,presence,of,both,modes,with,similar,damping,rates.,This,necessitates,in,some,cases,resolving,the,damping,of,both,modes,over,many,orders,of,magnitude,before,the,least,damped,mode,dominates,the,numerical,solution,–,a,challenging,task,given,the,ever-increasing,ﬁlamentation,of,the,velocity,space,due,to,phase,mixing.,This,should,be,possible,to,eliminate,by,carefully,initialising,the,simulation,so,that,only,the,least,damped,mode,is,present,,though,we,have,not,yet,attempted,this.,An,example,of,a,case,in,which,both,modes,are,present,,as,well,as,a,case,in,which,they,are,not,,is,given,in,Fig.,2.,00.511.522.500.511.52ω/kkvth(ni+nn)Rin/(cid:12)(cid:12)(cid:12)kkvth(cid:12)(cid:12)(cid:12)1/41/2124−1.4−1.2−1−0.8−0.6−0.4−0.200.511.52γ/kkvth(ni+nn)Rin/(cid:12)(cid:12)(cid:12)kkvth(cid:12)(cid:12)(cid:12)1/41/2124,Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,6,Figure,2.,Time,evolutions,of,the,absolute,value,of,the,spatially-varying,electrostatic,potential,δφ,(blue),and,the,result,of,a,least-squares,ﬁt,(orange),to,obtain,the,damping,rate,,frequency,and,phase.,The,left,plot,corresponds,to,Te/Ti,=,2,and,˜Rin,=,0,,and,the,right,plot,corresponds,to,Te/Ti,=,1,and,˜Rin,=,4,(,0.7,in,terms,of,the,normalised,frequency,used,in,Fig.,1).,≈,5.,Future,plans,With,the,core,1+1D,code,developed,and,successfully,benchmarked,against,analytical,theory,,we,will,now,turn,our,attention,to,including,kinetic,electron,dynamics,and,to,trialing,a,version,of,the,moment-based,approach,in,which,the,density,moment,is,split,In,parallel,,we,plan,to,explore,the,idea,of,using,an,oﬀ,from,the,kinetic,equation.,interpolation-free,semi-Lagrange,scheme,[10],for,the,solution,of,the,kinetic,equation,as,a,novel,approach,to,ameliorating,the,time,step,restriction,posed,by,the,kinetic,electron,dynamics.,Appendix,A.,Operator,splitting,We,evolve,Eqs.,(12)-(15),using,Heun’s,method,with,‘ﬂip-ﬂop’,Lie,operator,splitting.,To,facilitate,this,it,is,convenient,to,deﬁne,the,vector,f,whose,components,are,the,distribution,functions,˜fi,and,˜fn;,i.e.,,f,=,(,˜fi,,˜fn)T,.,In,terms,of,f,the,system,of,equations,is,∂f,∂t,+,A[f,],+,B[f,],=,C[f,],,where,and,A[f,],.,=,˜v,(cid:107),(cid:32),B[f,],.,=,−,1,2,(cid:32),1,0,0,vth,n/vth,i,(cid:33),∂,˜φ/∂,˜z,0,0,0,(cid:33),∂f,∂z,,,,,∂f,∂˜v,(cid:107),C[f,],.,=,˜Rin,−,(cid:32),˜ni,˜nn,−,˜ni,(cid:33),˜nn,−,f,.,(A.1),(A.2),(A.3),(A.4),Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,7,Splitting,the,operators,and,employing,Heun’s,method,results,in,the,following,time,advance,scheme:,f,(1),n,A,=,f,n,f,(2),n,A,=,f,n,1,2,f,n,A,=,∆tA,[f,n],−,∆tA,f,(1),n,A,−,(cid:104),A,+,f,(2),n,f,(1),n,A,A,−,(cid:16),f,(1),n,B,=,f,n,f,(2),n,B,=,f,n,A,−,1,2,f,n,B,=,(cid:16),∆tB,[f,n,A],f,(1),n,A,∆tB,(cid:104),B,+,f,(2),n,f,(1),n,B,,,(cid:105),(cid:17),(cid:105),,,(cid:17),and,B,−,f,(1),n,C,=,f,n,f,(2),n,C,=,f,n,B,−,1,2,f,n+1,=,∆tC,[f,n,B],f,(1),n,B,∆tC,(cid:105),(cid:104),f,(1),n,C,+,fC(2),,n,(cid:17),(cid:16),(A.5),(A.6),(A.7),,,with,the,order,of,operations,(A,,followed,by,B,,followed,by,C),reversed,every,time,step,to,obtain,second,order,acuracy,in,∆t.,Numerical,study,of,1D,drift,kinetic,models,with,periodic,boundary,conditions,8,[1],F.,I.,Parra,,M.,Barnes,,and,M.,R.,Hardman.,1d,drift,kinetic,models,with,periodic,boundary,conditions.,Excalibur/Neptune,Report,,1:2047357–TN–01–02,M1.1,,2021.,[2],C.-W.,Shu,and,S.,Osher.,Eﬃcient,implementation,of,essentially,non-oscillator,shock-capturing,schemes.,J.,Comp.,Phys.,,77:439–471,,1988.,[3],S.,Gottlieb,and,C.-W.,Shu.,Total,variation,diminishing,runge-kutta,methods.,Mathematics,of,Computation,,67:73–85,,1998.,[4],S.,Gottlieb,,C.-W.,Shu,,and,E.,Tadmor.,Strong,stability-preserving,high-order,time,discretization,methods.,SIAM,Rev.,,43:89,,2001.,[5],D.,R.,Durran.,Numerical,methods,for,ﬂuid,dynamics.,Springer,,2010.,[6],J.,P.,Boyd.,Chebyshev,and,Fourier,Spectral,Methods.,Dover,,New,York,,2000.,[7],M.,Abramowitz,and,I.,A.,Stegun.,Handbook,of,Mathematical,Functions,with,Formulas,,Graphs,,and,Mathematical,Tables.,Dover,,New,York,,1972.,[8],Matteo,Frigo,and,Steven,G.,Johnson.,The,design,and,implementation,of,FFTW3.,Proceedings,of,the,IEEE,,93(2):216–231,,2005.,Special,issue,on,“Program,Generation,,Optimization,,and,Platform,Adaptation”.,[9],C.,W.,Clenshaw,and,A.,R,Curtis.,A,method,for,numerical,integration,on,an,automatic,computer.,Numerische,Mathematik,,2:197,,1960.,[10],H.,Ritchie.,Eliminating,the,interpolation,associated,with,the,semi-lagrange,scheme.,Monthly,Weather,Review,,114:135–146,,1986.

:pdfembed:`src:_static/TN-02_Numerical1DDriftKineticPeriodicBoundaryConditions.pdf, height:1600, width:1100, align:middle`