TN-07_AHigherOrderFiniteElementImplementationFullFlandauFokkerPlanckCollisionOperatorC
======================================================================================

.. meta::
   :description: technical note

   :keywords: Report,2070839-TN-07,1,A,higher-order,finite-element,implementation,of,the,full-F,Landau,Fokker-Planck,collision,operator,for,charged,particle,collisions,in,a,low,density,plasma,M.,R.,Hardman1,2,,M.,Abazorius2,,J.,Omotani3,,M.,Barnes2,,S.,L.,Newton3,,and,F.,I.,Parra4,1,Tokamak,Energy,Ltd,,173,Brook,Drive,,Milton,Park,,Abingdon,,OX14,4SD,,United,Kingdom,2,Rudolf,Peierls,Centre,for,Theoretical,Physics,,University,of,Oxford,,Clarendon,Laboratory,,Parks,Road,,Oxford,OX1,3PU,,United,Kingdom,3,United,Kingdom,Atomic,Energy,Authority,,Culham,Science,Centre,,Abingdon,,Oxon,,OX14,3DB,,UK,4,Princeton,Plasma,Physics,Laboratory,,P.O.,Box,451,,Princeton,,New,Jersey,08540,,United,States,E-mail:,michael.hardman@tokamakenergy.co.uk,1.,Introduction,A,low,density,plasma,is,one,that,can,be,accurately,described,by,the,one-point,particle,probability,distribution,function,Fs(r,,v).,The,distribution,function,provides,us,with,the,probability,p,=,Fs(r,,v),d3rd3v,to,find,a,particle,of,species,s,in,the,phase,space,volume,around,the,phase,space,position,(r,,v),,with,r,the,particle,position,and,v,the,particle,velocity.,An,equation,for,the,distribution,function,may,be,obtained,from,the,BBGKY,hierarchy,[1,,2],,which,converts,an,N,-body,Hamiltonian,system,describing,a,plasma,or,gas,into,a,statistical,description.,The,resulting,equation,has,the,form,∂Fs,∂t,+,v,·,∇Fs,+,Zse,ms,(E,+,v,×,B),·,∂Fs,∂v,(cid:88),=,s′,Css′,[Fs,,Fs′],,,(1),where,the,left-hand,side,of,the,equation,is,the,Vlasov,equation,,accounting,for,the,acceleration,of,particles,by,the,large-scale,electromagnetic,fields,,and,the,Boltzmann,collision,operator,on,the,right-hand,side,of,the,equation,accounts,for,the,interactions,of,particles,of,species,s,with,local,small-scale,electromagnetic,fields,generated,by,interactions,between,particles,of,differing,species,s′,at,the,same,position,r.,Here,,ms,is,the,species,mass,,Zs,is,the,species,charge,number,,e,is,the,unit,charge,,and,E,and,B,are,the,electric,and,magnetic,fields,,respectively,,If,the,interaction,cross,section,is,chosen,to,be,the,1/r2,electrostatic,potential,,then,the,collision,operator,becomes,the,well-known,Landau,Fokker-Planck,collision,operator,Report,2070839-TN-07,[3,,4,,5,,6]:,Css′,[Fs,,Fs′],=,(cid:26)(cid:90),γss′,ms,∂,∂v,·,∂2g,∂v∂v,(cid:20),Fs′(v′),ms,·,∂Fs,∂v,−,Fs(v),ms′,∂Fs′,∂v′,(cid:21),(cid:27),,,d3v′,where,γss′,=,2πZ,2,2e4,ln,Λss′,s,Zs′,(4πϵ0)2,,,with,ln,Λss′,the,Coulomb,logarithm,[3,,6,,4,,5],,and,g,=,|v,−,v′|.,2,(2),(3),(4),The,operator,is,integro-differential,because,the,effect,of,small-angle,collisions,dominate,over,the,rarer,large-angle,collisions.,The,complex,structure,of,the,collision,operator,obscures,four,key,properties,that,we,note.,First,,the,collision,operator,conserves,particle,density,,i.e.,,(cid:90),Css′,[Fs,,Fs′],d3v,=,0.,Second,,the,collision,operator,conserves,the,total,momentum,in,a,collision,,i.e.,,(cid:90),(msv,Css′,[Fs,,Fs′],+,ms′v,Cs′s,[Fs′,,Fs]),d3v,=,0.,The,same,is,true,for,the,total,energy:,(cid:90),(cid:18),1,2,ms|v|2,Css′,[Fs,,Fs′],+,ms′|v|2,Cs′s,[Fs′,,Fs],(cid:19),d3v,=,0.,1,2,(5),(6),(7),Finally,,Boltzmann’s,H-theorem,applied,to,same-species,collisions,[5,,6,,2],proves,that,the,entropy,production,(cid:90),˙Ss,=,−,ln,Fs,Css,[Fs,,Fs],d3v,≥,0,(8),with,equality,if,and,only,if,Fs,is,a,Maxwellian,distribution,described,by,the,local,density,ns,,mean,velocity,us,,and,temperature,Ts,,i.e.,,Fs,=,F,M,s,=,ns,π3/2v2,th,s,(cid:34),exp,−,(cid:18),v,−,us,vth,s,(cid:19)2(cid:35),,,(9),with,vth,s,=,(cid:112)2Ts/ms.,Implementing,the,full,Landau,collision,operator,is,challenging,because,of,the,nonlinear,and,integro-differential,nature,of,the,operator.,For,a,given,distribution,function,Fs,,we,must,carry,out,a,series,of,difficult,integrals,to,find,the,coefficients,of,the,operator.,Whilst,previous,authors,have,implemented,the,full,operator,,see,,e.g.,,[7,,8,,9,,10],,including,implementations,of,the,underlying,Boltzmann,operator,[11],,it,is,more,typical,to,use,asymptotic,expansions,to,linearised,the,operator,for,its,use,in,a,Report,2070839-TN-07,3,specific,application,(e.g.,transport,theory,or,collisional,closures,[4,,5,,12,,13]).,It,is,also,common,to,write,down,an,ad-hoc,diffusive,model,operator,which,may,be,solved,rapidly,yet,still,has,the,conservation,or,H-theorem,properties,desired,for,the,physics,of,interest,[14,,15,,16,,17].,In,core-transport,plasma,turbulence,applications,,collisional,relaxation,timescales,are,typically,long,compared,to,the,non-linear,turnover,time,of,the,turbulent,eddies,,meaning,that,the,linearised,Landau,operator,or,ad-hoc,model,operators,are,an,appropriate,cost,saving,device:,physically,,all,these,operators,are,required,to,do,is,dissipate,fine,velocity-space,structure,[15,,16].,The,system,is,known,to,be,approximately,in,thermal,equilibrium,because,the,system,is,closed,[13,,18],,meaning,that,the,distribution,function,is,never,far,from,the,Maxwellian,around,which,the,collision,operator,is,easily,linearised.,However,,in,the,scrape-off,layer,the,distribution,function,of,the,bulk,plasma,may,well,be,far,from,Maxwellian.,This,is,due,to,the,presence,of,the,divertor,plate,or,limiter,,which,makes,the,system,open,,preventing,local,thermal,equilibrium.,In,addition,,hot,particles,may,transit,rapidly,from,the,core,to,the,open,field,lines,at,the,edge,where,the,bulk,of,the,plasma,is,expected,to,be,cooler,,potentially,resulting,in,a,bimodal,distribution,of,particle,energies.,In,general,,the,exact,steady,state,distribution,is,not,known.,In,summary,,it,is,not,clear,whether,or,not,a,model,or,linearised,collision,operator,is,adequate,for,modelling,the,plasma,on,open,field,lines.,The,only,rigorous,choice,of,operator,is,the,full-F,Landau,operator.,In,this,report,we,describe,the,implementation,of,the,full-F,Landau,operator,using,higher-order,finite,elements,for,gyrotropic,distribution,functions.,Gyrotropic,distributions,are,independent,of,the,gyrophase,angle,ϑ,that,measures,the,position,of,the,particle,in,the,plane,perpendicular,to,the,magnetic,field.,Mathematically,,this,means,that,we,support,Fs,=,Fs(v∥,,v⊥),,with,the,cylindrical,velocity,space,coordinates,(v∥,,v⊥,,ϑ),defined,by,v∥,=,v,·,b,,v⊥,=,|v,−,v∥b|,,tan,ϑ,=,−,v,·,e2,v,·,e1,,,or,equivalently,,v,=,v∥b,+,v⊥(cos,ϑe1,−,sin,ϑe2).,(10),(11),The,basis,vector,b,=,B/|B|,is,the,unit,vector,in,the,direction,of,the,magnetic,field.,The,vectors,e1,and,e2,are,orthogonal,to,b,and,satisfy,b,·,e1,×,e2,=,1,,e1,·,b,=,0,,e2,·,b,=,0.,(12),The,numerical,implementation,described,in,this,report,ensures,the,satisfaction,of,the,conservative,properties,(5)-(7),by,achieving,high,accuracy,with,the,weak,formulations,and,adequate,numerical,resolution.,Almost,exact,numerical,conservation,is,obtained,with,ah-hoc,numerical,conserving,terms,,which,are,necessary,to,preserve,a,stable,solution,at,long,times.,To,avoid,carrying,out,costly,numerical,integration,to,find,the,coefficients,in,the,whole,of,the,velocity,space,,we,use,the,Rosenbluth-MacDonald-,Judd,(RMJ),form,of,the,collision,operator,(which,is,a,Fokker-Planck,form),[3],and,we,Report,2070839-TN-07,4,solve,elliptic,PDEs,for,the,required,coefficients,,called,Rosenbluth,potentials,,using,the,higher-order,finite,element,method,and,boundary,conditions,obtained,from,the,Green’s,functions,at,the,limits,of,the,velocity,space.,This,numerical,strategy,optimises,the,scheme,for,scalability.,The,remainder,of,the,report,is,structured,as,follows.,In,the,next,section,,we,write,the,collision,operator,in,the,RMJ,form.,In,section,3,,we,obtain,the,weak-form,representation,of,the,problem,that,we,will,implement,numerically.,In,section,4,we,discuss,the,boundary,conditions,imposed,on,Fs,during,an,explicit,time,advance.,In,section,5,we,prescribe,numerical,conserving,terms,for,use,in,the,time,advance.,In,section,6,we,provide,numerical,results,from,testing,our,implementation.,In,section,7,,we,discuss,the,outlook,for,the,use,of,the,operator,in,a,production,code.,Finally,,Appendix,A,and,Appendix,B,contain,useful,results,pertaining,to,numerical,integration,of,the,Green’s,functions,of,the,derivatives,of,the,Rosenbluth,potentials.,2.,Rosenbluth-MacDonald-Judd,form,of,the,collision,operator,The,operator,in,the,Rosenbluth-MacDonald-Judd,(RMJ),form,[3,,6],coordinates,is,most,usefully,written,in,terms,of,collisional,fluxes:,in,(v∥,,v⊥),Css′,[Fs,,Fs′],=,∂Γ∥,∂v∥,+,1,v⊥,∂,∂v⊥,(v⊥Γ⊥),.,where,the,fluxes,are,defined,by,and,Γ∥,=,γss′,m2,s,(cid:18),∂Fs,∂v∥,∂2Gs′,∂v∥,2,+,∂Fs,∂v⊥,∂2Gs′,∂v⊥∂v∥,−,2,ms,ms′,Fs,∂Hs′,∂v∥,Γ⊥,=,γss′,m2,s,(cid:18),∂Fs,∂v∥,∂2Gs′,∂v∥∂v⊥,+,∂Fs,∂v⊥,∂2Gs′,∂v⊥,2,−,2,ms,ms′,Fs,∂Hs′,∂v⊥,(cid:19),,,(cid:19),,,and,the,Rosenbluth,potentials,are,and,(cid:90),Gs′(v),=,Fs′(v′)g,d3v′,Hs′(v),=,(cid:90),Fs′(v′),g,d3v′.,(13),(14),(15),(16),(17),In,the,drift-kinetic,limit,the,largest,piece,of,the,distribution,functions,are,independent,of,gyroangle,[4,,19],,i.e.,,Fs,=,Fs(v∥,,v⊥),and,Fs′,=,Fs′(v∥,,v⊥).,In,terms,of,(v∥,,v⊥),coordinates,,for,gyrotropic,distributions,the,Rosebluth,potentials,simplify,to,Gs′,=,(cid:90),∞,(cid:90),∞,0,−∞,(cid:16)(cid:0)v∥,−,v′,∥,4,(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)1/2,E(m(v∥,,v⊥,,v′,∥,,v′,⊥))Fs′(v′,∥,,v′,⊥)v′,⊥,dv′,∥dv′,⊥,,(18),Report,2070839-TN-07,5,and,Hs′,=,where,(cid:90),∞,(cid:90),∞,0,−∞,(cid:16)(cid:0)v∥,−,v′,∥,4,(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)−1/2,K(m(v∥,,v⊥,,v′,∥,,v′,⊥))Fs′(v′,∥,,v′,⊥)v′,⊥,dv′,∥dv′,⊥,,(19),(20),(21),(22),(23),(24),m(v∥,,v⊥,,v′,∥,,v′,⊥),=,4v⊥v′,⊥,(cid:16)(cid:0)v∥,−,v′,∥,(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)−1,,,and,we,have,used,the,definitions,of,the,complete,elliptic,integral,of,the,first,kind,K(m),=,(cid:90),π/2,0,1,1,−,m,sin2,θ,(cid:112),dθ,and,the,complete,elliptic,integral,of,the,second,kind,E(m),=,(cid:90),π/2,(cid:112),0,1,−,m,sin2,θ,dθ.,2.1.,Finding,elliptic,problems,for,the,Rosenbluth,potentials,As,noted,in,the,original,derivation,by,Rosenbluth,,MacDonald,,and,Judd,,[3],,the,potentials,defined,by,eqations,(16),and,(17),may,also,be,defined,as,the,solutions,of,the,elliptic,problems,∂2G,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:19),(cid:18),v⊥,∂G,∂v⊥,=,2H,,and,∂2H,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:19),(cid:18),v⊥,∂H,∂v⊥,=,−4πF.,We,recognise,the,Poisson’s,equation,for,the,Rosenbluth,potential,H,,and,the,biharmonic,equation,for,G.,Obtaining,the,Rosenbluth,potentials,through,an,elliptic,solve,is,potentially,numerically,advantageous,compared,to,evaluating,the,Green’s,function,directly.,This,is,because,sparse,matrix,methods,can,be,used,to,invert,the,Laplacian,operator,,giving,an,operation,scaling,of,N,2,,where,N,is,the,size,of,one,of,the,velocity,dimensions,,whereas,a,direct,evaluation,of,the,Green’s,function,leads,to,a,scaling,of,the,order,of,N,4.,An,important,detail,is,that,for,a,finite,simulation,domain,,adequate,boundary,conditions,must,be,supplied,when,carrying,out,the,elliptic,solve.,In,practice,this,necessitates,the,order,N,3,operation,of,obtaining,the,boundary,data,through,the,Green’s,function.,An,order,N,improvement,in,the,overall,scaling,is,not,to,be,dismissed,,and,parallelisation,over,many,processes,may,be,able,to,alleviate,a,further,order,N,(corresponding,to,the,size,of,the,boundary,in,velocity,space),as,evaluating,the,boundary,data,is,embarrassingly,parallel.,Having,motivated,finding,the,coefficients,appearing,in,the,collision,operator,through,elliptic,solves,,it,remains,to,find,the,appropriate,PDEs.,Direct,differentiation,of,equations,(23),and,(24),yields,the,required,differential,definitions,of,the,coefficients,,as,we,now,demonstrate.,We,can,solve,for,G20,=,∂2G,2,∂v∥,(25),Report,2070839-TN-07,by,simply,differentiating,equation,(23),to,find,∂2G20,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:19),(cid:18),v⊥,∂G20,∂v⊥,=,2,∂2H,2,.,∂v∥,The,next,coefficient,that,we,require,is,G01,=,∂G,∂v⊥,.,6,(26),(27),To,obtain,an,elliptic,equation,for,G01,we,differentiate,equation,(23),with,respect,to,v⊥,to,find,that,∂2G01,∂v∥,2,+,∂2G01,∂v⊥,2,+,∂,∂v⊥,(cid:19),(cid:18),G01,v⊥,=,2,∂H,∂v⊥,.,Using,the,identity,∂2G,∂v⊥,2,=,1,v⊥,∂,∂v⊥,(cid:19),(cid:18),v⊥,∂G,∂v⊥,−,1,v⊥,∂G,∂v⊥,,,we,can,rewrite,the,equation,for,G01,in,a,form,amenable,to,integration,by,parts:,∂2G01,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:18),v⊥,∂G01,∂v⊥,(cid:19),−,G01,v2,⊥,=,2,∂H,∂v⊥,.,The,next,coefficient,that,we,require,is,G11,=,∂2G,∂v∥∂v⊥,.,(28),(29),(30),(31),The,differential,equation,that,defines,G11,may,be,obtained,from,(30),by,direct,differentiation.,The,result,is,∂2G11,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:18),v⊥,∂G11,∂v⊥,(cid:19),−,G11,v2,⊥,=,2,∂2H,∂v∥∂v⊥,.,To,complete,the,set,of,coefficients,derived,from,G,,an,algebraic,equation,for,G02,=,∂2G,2,∂v⊥,(32),(33),can,be,obtained,from,the,defining,equation,for,G,,once,the,other,derivatives,are,obtained.,We,have,that,∂2G,∂v⊥,2,=,2H,−,∂2G,∂v∥,2,−,1,v⊥,∂G,∂v⊥,or,in,the,notation,introduced,here,G02,=,2H,−,G20,−,G01,v⊥,.,,,(34),(35),Report,2070839-TN-07,7,It,is,also,possible,to,obtain,a,convenient,elliptic,problem,for,G02,by,differentiation,of,(23),with,appropriate,use,of,the,identity,(35).,We,implement,∂2G02,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:18),v⊥,∂G02,∂v⊥,(cid:19),−,4G02,v2,⊥,=,2,=,2,−,(cid:18),∂2H,∂v⊥,∂,2,∂v⊥,v⊥,4H,v2,⊥,+,v⊥,∂H,∂v⊥,2,v2,⊥,(cid:19),∂2G,2,∂v∥,2,v⊥,−,∂H,∂v⊥,−,4H,v2,⊥,+,2,v2,⊥,∂2G,2,∂v∥,(36),Note,that,solving,for,G02,requires,knowledge,of,H,and,other,derivatives,of,G,,regardless,of,the,formulation.,The,primary,benefit,of,using,the,solution,of,the,elliptic,equation,(30),rather,than,the,algebraic,equation,(35),is,numerical,accuracy.,In,analogy,to,the,derivations,above,,it,is,straightforward,to,show,that,the,equations,for,are,given,by,H10,=,∂H,∂v∥,,,and,H01,=,∂H,∂v⊥,∂2H10,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:19),(cid:18),v⊥,∂H10,∂v⊥,=,−4π,∂F,∂v∥,.,(37),(38),and,(cid:18),(cid:19),−,v⊥,2,+,1,v⊥,∂2H01,∂v∥,=,−4π,∂,∂v⊥,H01,v2,⊥,∂H01,∂v⊥,respectively.,Note,that,we,have,written,the,elliptic,equations,(26),,(30),,(32),,(36),,(38),,and,(39),,in,a,form,that,will,be,amenable,to,integration,by,parts,in,the,following,test-function,analysis,required,for,a,weak-form,implementation.,The,numerical,implementation,of,these,equations,first,uses,equations,(24),,(38),and,(39),to,find,H,and,its,derivatives,from,F,.,Then,,equations,(26),,(30),,(32),,and,(36),may,be,solved,for,the,derivatives,of,G.,∂F,∂v⊥,(39),,,2.2.,Evaluating,the,boundary,data,Numerical,integration,of,diverging,functions,is,challenging.,To,obtain,the,boundary,data,required,to,solve,the,elliptic,problems,(26),,(30),,(32),,(38),,and,(39),,we,aim,to,avoid,Green’s,functions,with,denominators,that,have,large,powers,of,g.,We,can,achieve,this,by,using,integration,by,parts,to,relax,the,divergence.,This,has,the,effect,of,making,the,required,Green’s,functions,similar,in,form,to,those,required,for,G,and,H.,We,start,by,computing,(cid:90),=,∂Gs′,∂v,Fs′(v′),∂g,∂v,(cid:90),d3v′,=,−,Fs′(v′),∂g,∂v′,d3v′,,(40),where,we,have,used,that,∂g/∂v,=,−∂g/∂v′.,Using,integration,by,parts,,and,that,Fs′(v′),→,0,as,|v′|,→,0,,we,find,that,∂Gs′,∂v,=,(cid:90),∂Fs′,∂v′,g,d3v′.,(41),Report,2070839-TN-07,We,can,carry,out,this,procedure,a,second,time,to,find,that,∂2Gs′,∂v∂v,=,(cid:90),∂2Fs′,∂v′∂v′,g,d3v′.,We,can,use,the,same,method,to,find,that,∂Hs′,∂v,=,(cid:90),∂Fs′,∂v′,1,g,d3v′.,Without,the,second,integration,by,parts,,we,also,have,∂2Gs′,∂v∂v,=,(cid:90),∂Fs′,∂v′,v,−,v′,g,d3v′.,8,(42),(43),(44),Equations,(41),and,(43),are,vector,equation,and,(42),and,(44),are,tensor,equations.,We,extract,the,relevant,coefficients,by,taking,the,dot,product,with,the,unit,vectors,b,and,e⊥,,noting,that,Assuming,that,Fs′,=,Fs′(v′,⊥),(which,implies,that,Gs′,=,Gs′(v∥,,v⊥),,we,find,that,and,∂,∂v,=,b,+,e⊥,∂,∂v∥,∂,∂v⊥,e⊥,×,b,v⊥,∂,∂ϑ,,,∂,∂v′,∥,+,e′,⊥,∂,∂v′,⊥,e′,⊥,×,b,v′,⊥,∂,∂ϑ′,.,+,+,∂,∂v′,=,b,∥,,v′,(cid:90),∂Fs′,∂v′,⊥,(e⊥,·,e′,⊥)g,d3v′,=,2π,(cid:90),(cid:90),∂Fs′,∂v′,⊥,IG1,v′,⊥dv′,⊥dv′,∥,,d3v′,=,2π,IH0,v′,⊥dv′,⊥dv′,∥,,=,be⊥,:,∂2Gs′,∂v∂v,(cid:90),=,∂2Fs′,∥∂v′,∂v′,⊥,(e⊥,·,e′,⊥)g,d3v′,(cid:90),(cid:90),=,2π,∂2Fs′,∥∂v′,∂v′,⊥,1,g,(cid:90),∂Fs′,∂v′,∥,e⊥,·,e′,⊥,g,(cid:90),∂Fs′,∂v′,⊥,∂2Gs′,∂v∂v,(cid:90),(cid:90),∂Fs′,∂v′,∥,=,=,2π,∂2Gs′,∂v∥,2,=,bb,:,IG1,v′,⊥dv′,⊥dv′,∥,,(cid:90),(cid:90),∂Fs′,∂v′,∥,(cid:90),(cid:90),∂Fs′,∂v′,⊥,(cid:90),∂Fs′,∂v′,∥,v∥,−,v′,∥,g,d3v′,(v∥,−,v′,∥)IH0,v′,⊥dv′,⊥dv′,∥,,d3v′,=,2π,IH1,v′,⊥dv′,⊥dv′,∥.,∂Gs′,∂v⊥,=,e⊥,·,∂Gs′,∂v,=,∂2Gs′,∂v∥∂v⊥,∂Hs′,∂v∥,=,b,·,∂Hs′,∂v,=,∂Hs′,∂v⊥,=,e⊥,·,∂Hs′,∂v,=,and,∂2Gs′,∂v⊥,2,=,e⊥e⊥,:,∂2Gs′,∂v∂v,(cid:90),(cid:90),∂Fs′,∂v′,⊥,=,2π,=,(cid:90),∂Fs′,∂v′,⊥,(e⊥,·,e′,⊥),v⊥,−,v′,⊥(e⊥,·,e′,g,⊥),d3v′,(v⊥IH1,−,v′,⊥IH2),v′,⊥dv′,⊥dv′,∥,,(45),(46),(47),(48),(49),(50),(51),(52),9,(53),(54),(55),(56),Report,2070839-TN-07,where,IG1,=,1,2π,(cid:90),π,g,(e⊥,·,e′,⊥),dϑ′,,−π,1,2π,(cid:90),π,−π,(cid:90),π,dϑ′,,1,g,−π,e⊥,·,e′,⊥,g,dϑ′,,IH0,=,IH1,=,1,2π,and,IH2,=,(cid:90),π,1,2π,⊥)2,(e⊥,·,e′,g,dϑ′.,−π,formulation,is,the,integrands,only,diverge,The,main,advantage,of,logarithmically,where,v′,⊥,=,v⊥.,This,kind,of,divergence,can,be,handled,numerically,by,a,change,of,variables,in,the,affected,elements,[20].,The,functions,IG1,,IH0,,IH1,,and,IH2,are,evaluated,in,Appendix,A.,this,∥,=,v∥,and,v′,that,3.,Obtaining,the,weak,formulation,of,the,problem,We,consider,the,problem,∂F,∂t,=,∂Γ∥,∂v∥,+,1,v⊥,∂,∂v⊥,(v⊥Γ⊥),.,(57),in,v∥,∈,[−L∥,,L∥],,v⊥,∈,[0,,L⊥],and,t,,where,L∥,and,L⊥,are,the,maximum,values,of,v∥,and,v⊥,on,the,grid,,respectively.,The,solution,F,=,F,(v∥,,v⊥,,t),,and,the,coefficients,Γ∥,=,Γ∥(v∥,,v⊥,,t),=,Γ∥[F,(v∥,,v⊥,,t)],and,Γ⊥,=,Γ⊥(v∥,,v⊥,,t),=,Γ⊥[F,(v∥,,v⊥,,t)],are,functionals,of,F,.,We,note,that,the,fluxes,in,velocity,space,are,defined,explicitly,by,equations,(14),and,(15).,We,divide,the,domain,into,a,rectangular,grid,of,N2D,=,Nelement,∥Nelement,⊥,elements.,We,use,Nelement,∥,1D,elements,in,the,v∥,direction,and,Nelement,⊥,1D,elements,in,the,v⊥,direction.,Each,2D,element,is,an,outer,product,of,two,1D,elements.,On,each,1D,element,we,represent,the,function,with,Lagrange,polynomials,of,order,Ngrid,+,1,using,the,(normalised),points,within,the,elements,xj,∈,{x0,,x1,,...,,xNgrid−1,,xNgrid},(58),with,x0,=,−1,and,xNgrid,=,1,(Lobatto,points),on,elements,that,do,not,include,v⊥,=,0.,On,the,element,including,v⊥,=,0,,we,take,xNgrid,=,1,but,we,use,Radau,quadrature,to,ensure,that,x0,>,−1.,The,transformation,between,(v∥,,v⊥),and,the,local,coordinate,x(r),in,the,rth,1D,element,is,v∥,=,s(r),∥,x(r),+,c(r),∥,,,v⊥,=,s(r),⊥,x(r),+,c(r),⊥,(59),∥,,,c(r),∥,,,s(r),⊥,and,c(r),where,s(r),⊥,are,constants,in,each,element,(labelled,here,by,r),which,may,vary,between,elements,,and,x(r),∈,[−1,,1],for,all,r,(saving,the,element,including,the,origin,of,v⊥,,which,has,x(r),∈,(−1,,1]).,Report,2070839-TN-07,3.1.,The,basis,functions,We,introduce,2D,basis,functions,Φ(rp),jk,(cid:0)v∥,,v⊥,(cid:1),=,φ(r),j,(cid:0)v∥,(cid:1),φ(p),k,(v⊥),,,where,the,1D,basis,functions,are,φ(r),j,(v),=,lj,(cid:0)x(r),(v)(cid:1),H,(cid:16),v,−,v,(cid:16),x(r),0,(cid:17)(cid:17),(cid:16),(cid:16),v,H,x(r),Ngrid,(cid:17),(cid:17),,,−,v,10,(60),(61),where,v,is,a,placeholder,for,either,v∥,or,v⊥,and,lj,is,the,jth,Lagrange,polynomial,on,the,element.,Expanding,the,solution,in,these,basis,functions,,we,write,F,(v∥,,v⊥),=,=,(cid:88),(cid:88),r,p,(cid:88),j,k,(cid:88),r,p,j,k,jk,Φ(rp),F,rp,jk,(v∥,,v⊥),jk,φ(r),F,rp,j,(v∥)φ(p),k,(v⊥),,with,F,rp,jk,=,F,(cid:16),(cid:16),v∥,x(r),j,(cid:17),,,v⊥,(cid:17)(cid:17),(cid:16),x(p),k,.,Note,that,the,basis,functions,have,the,cardinality,property,φ(r),j,(cid:16),v,(cid:16),x(p),k,(cid:17)(cid:17),=,δjkδrp.,3.2.,The,projection,(62),(63),(64),To,project,equation,(57),onto,the,basis,functions,Φ(rp),jk,basis,function,Φ(qs),mn,2D,element.,The,limits,of,this,element,are,v(q),(cid:1),,we,multiply,by,the,(cid:1),,and,integrate,over,velocity,space,corresponding,to,a,single,(cid:17),,,(cid:0)v∥,,v⊥,(cid:0)v∥,,v⊥,,,v(q),(cid:17),(cid:16),(cid:16),∥l,=,v∥,x(q),0,∥u,=,v∥,x(q),Ngrid,v(s),⊥u,=,v⊥,(cid:16),x(s),Ngrid,(cid:17),,,and,v(s),⊥l,=,v⊥,(cid:17),(cid:16),x(s),0,,,respectively.,Then,we,have,that,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,Φ(qs),mn,∂F,∂t,v⊥dv⊥dv∥,=,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,Φ(qs),mn,(cid:18),∂Γ∥,∂v∥,+,1,v⊥,∂,∂v⊥,(cid:19),(v⊥Γ⊥),v⊥dv⊥dv∥.,(65),Report,2070839-TN-07,3.3.,The,mass,matrix,The,left,hand,side,of,equation,(65),takes,the,form,11,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,Φ(qs),mn,∂F,∂t,=,(cid:88),(cid:88),∂F,rp,jk,∂t,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,mn,Φ(rp),Φ(qs),jk,v⊥dv⊥dv∥,j,k,∂F,qs,jk,∂t,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,mn,Φ(qs),Φ(qs),jk,v⊥dv⊥dv∥,∂F,qs,jk,∂t,(cid:90),v(q),∥u,v(q),∥l,m,(v∥)φ(q),φ(q),j,(v∥),dv∥,(cid:90),v(s),⊥u,v(s),⊥l,n,(v⊥)φ(s),φ(s),k,(v⊥),v⊥dv⊥,M,(q),∥mjM,(s),⊥nk,∂F,qs,jk,∂t,.,r,p,(cid:88),=,j,k,(cid:88),j,k,(cid:88),j,k,=,=,We,deal,with,the,separated,integrals,as,follows,M,(q),∥mj,=,(cid:90),v(q),∥u,v(q),∥l,m,(v∥)φ(q),φ(q),j,(v∥),dv∥,=,s(q),∥,(66),(cid:90),1,−1,lm,(x),lj,(x),dx,(67),M,(s),⊥nk,=,(cid:90),v(s),⊥u,v(s),⊥l,n,(v⊥)φ(s),φ(s),k,(v⊥),v⊥dv⊥,=,s(s),⊥,(cid:90),1,−1,ln,(x),lk,(x),(cid:16),⊥,x,+,c(s),s(s),⊥,(cid:17),dx,(68),The,product,M,(q),∥mjM,(s),⊥nk,acting,on,F,rp,jk,is,often,referred,to,as,a,mass,matrix.,3.4.,The,nonlinear,stiffness,matrices,for,the,collision,operator,The,form,of,the,right,hand,side,of,equation,(65),and,the,forms,of,the,fluxes,,given,by,equations,(14),and,(15),,respectively,,suggest,that,we,should,integrate,by,parts,to,bring,all,derivatives,down,to,first,order.,Carrying,out,this,step,,we,find,that,for,the,parallel,flux,term,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,Φ(qs),mn,∂Γ∥,∂v∥,v⊥dv⊥dv∥,=,(cid:34)(cid:90),v(q),⊥u,v(q),⊥l,(cid:35)v(s),∥u,Φ(qs),mn,Γ∥,v⊥dv⊥,−,v(s),∥l,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,∂Φ(qs),mn,∂v∥,Γ∥,v⊥dv⊥dv∥,=,δmNgrid,(cid:90),v(q),⊥u,v(q),⊥l,mn,(v(s),Φ(qs),∥u,,,v⊥)Γ∥(v(s),∥u,,,v⊥),v⊥dv⊥,−,δm0,mn,(v(s),Φ(qs),∥l,,,v⊥)Γ∥(v(s),∥l,,,v⊥),v⊥dv⊥,(cid:90),v(q),⊥u,v(q),⊥l,(cid:90),v(s),⊥u,(cid:90),v(q),∥u,v(q),∥l,v(s),⊥l,−,∂Φ(qs),mn,∂v∥,Γ∥,v⊥dv⊥dv∥.,(69),Report,2070839-TN-07,Similarly,,for,the,perpendicular,flux,term,,we,have,that,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,Φ(qs),mn,v⊥,∂,∂v⊥,(v⊥Γ⊥),v⊥dv⊥dv∥,=,=,(cid:90),v(q),∥u,v(q),∥l,(cid:2)Φ(qs),mn,v⊥Γ⊥,(cid:3)v(s),⊥u,v(s),⊥l,dv∥,−,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,∂Φ(qs),mn,∂v⊥,Γ⊥,v⊥dv⊥dv∥,=,δnNgrid,(cid:90),v(q),∥u,v(q),∥l,mn,(v∥,,v(s),Φ(qs),⊥u)v(s),⊥uΓ⊥(v∥,,v(s),⊥u),dv∥,12,(70),−,δn0,(cid:90),v(q),∥u,v(q),∥l,mn,(v∥,,v(s),Φ(qs),⊥l,)v(s),⊥l,Γ⊥(v∥,,v(s),⊥l,),dv∥,−,(cid:90),v(q),∥u,(cid:90),v(s),⊥u,v(q),∥l,v(s),⊥l,∂Φ(qs),mn,∂v⊥,Γ⊥,v⊥dv⊥dv∥.,The,boundary,flux,terms,in,equations,(69),and,(70),will,cancel,identically,at,the,assembly,stage,,or,vanish,at,v⊥,=,0,or,v⊥,=,L⊥,,v∥,=,−L∥,,and,v∥,=,L∥,by,the,boundary,conditions,that,F,→,0,as,|v|,→,∞.,We,are,now,in,a,position,to,write,down,the,final,matrix,row.,We,use,the,expansion,(62),for,both,the,distribution,function,F,and,the,coefficients,that,are,derived,from,the,Rosenbluth,potentials,Gs′,and,Hs′.,With,these,choices,,recalling,the,definitions,of,the,fluxes,Γ∥,and,Γ⊥,,equations,(14),and,(15),,respectively,,the,result,is,M,(q),∥mjM,(s),⊥nk,(cid:88),j,k,∂F,qs,jk,∂t,=,−,γss′,m2,s,(cid:88),jklr,F,qs,jk,+,+,(cid:20),∂2Gs′,∂v⊥∂v∥,(cid:20),∂2Gs′,∂v⊥∂v∥,(cid:21)(qs),lr,(cid:21)(qs),lr,(cid:32)(cid:20)∂2Gs′,∂v∥,2,(cid:21)(qs),lr,∥2mjlY,(s),Y,(q),⊥0nkr,∥1mjlY,(s),Y,(q),⊥3nkr,−,2,∥3mjlY,(s),Y,(q),⊥1nkr,+,(cid:20),∂Hs′,∂v∥,(cid:21)(qs),ms,ms′,(cid:20),∂2Gs′,∂v⊥,2,(cid:21)(qs),lr,∥1mljY,(s),Y,(q),⊥0nkr,∥0mjlY,(s),Y,(q),⊥2nkr,lr,(71),−2,ms,ms′,(cid:20)∂Hs′,∂v⊥,(cid:21)(qs),lr,(cid:33),∥0mjlY,(s),Y,(q),⊥1nkr,+,δnNgrid,(cid:90),v(q),∥u,v(q),∥l,mn,(v∥,,v(s),Φ(qs),⊥u)v(s),⊥uΓ⊥(v∥,,v(s),⊥u),dv∥,−,δn0,(cid:90),v(q),∥u,v(q),∥l,mn,(v∥,,v(s),Φ(qs),⊥l,)v(s),⊥l,Γ⊥(v∥,,v(s),⊥l,),dv∥,+,δmNgrid,(cid:90),v(s),⊥u,v(s),⊥l,mn,(v(q),Φ(qs),∥u,,,v⊥)Γ∥(v(q),∥u,,,v⊥),v⊥dv⊥,−,δm0,(cid:90),v(s),⊥u,v(s),⊥l,mn,(v(q),Φ(qs),∥l,,,v⊥)Γ∥(v(q),∥l,,,v⊥),v⊥dv⊥,,Report,2070839-TN-07,where,we,have,defined,stiffness,matrices,with,three,indices,Y,(q),∥0mjl,=,Y,(q),∥2mjl,=,(cid:90),v(q),∥u,v(q),∥l,(cid:90),v(q),∥u,v(q),∥l,m,φ(q),φ(q),j,φ(q),l,dv∥,,Y,(q),∥1mjl,=,(cid:90),v(q),∥u,v(q),∥l,∂φ(q),m,∂v∥,j,φ(q),φ(q),l,dv∥,,∂φ(q),m,∂v∥,∂φ(q),j,∂v∥,φ(q),l,dv∥,,Y,(q),∥3mjl,=,(cid:90),v(q),∥u,v(q),∥l,φ(q),m,∂φ(q),j,∂v∥,φ(q),l,dv∥,,and,Y,(s),⊥0nkr,=,Y,(s),⊥2nkr,=,(cid:90),v(s),⊥u,v(s),⊥l,(cid:90),v(s),⊥u,v(s),⊥l,n,φ(q),φ(q),k,φ(q),r,v⊥,dv⊥,,Y,(s),⊥1nkr,=,(cid:90),v(s),⊥u,v(s),⊥l,∂φ(q),n,∂v⊥,φ(q),k,φ(q),r,v⊥,dv⊥,,∂φ(q),n,∂v⊥,∂φ(q),k,∂v⊥,φ(q),r,v⊥,dv⊥,,Y,(s),⊥3nkr,=,(cid:90),v(s),⊥u,v(s),⊥l,φ(q),n,∂φ(q),k,∂v⊥,φ(q),r,v⊥,dv⊥.,13,(72),(73),Note,that,the,stiffness,matrices,in,(71),are,all,1D,integrals,of,1D,basis,functions,,as,a,result,of,the,choice,to,use,the,representation,(62),where,the,2D,basis,function,Φ(qs),mn,(v∥,,v⊥),is,a,product,of,two,1D,lagrange,polynomials,–,one,for,the,v∥,dimension,,and,one,for,the,v⊥,dimension.,The,assembly,step,is,carried,out,by,defining,a,compound,index,that,indexes,over,all,the,rows,in,(v∥,,v⊥),,and,then,forming,a,matrix,equation,in,that,index.,We,use,continuity,of,F,to,demand,that,F,q,s,,,and,remove,the,duplicated,points,at,interior,element,boundaries,by,summing,the,matrix,rows,there.,Ngridk,=,F,q+1,s,=,F,q,s+1,j0,,,F,q,s,jNgrid,0k,3.5.,The,weak,form,of,the,equations,for,the,Rosenbluth,potentials,We,still,need,to,determine,the,coefficients,derived,from,the,Rosenbluth,potentials.,We,start,by,considering,the,solution,of,Poisson’s,equation,,equation,(24).,Multiplying,by,the,2D,basis,function,Φ(rp),k,(v⊥),and,integrating,over,velocity,space,we,have,jk,=,φ(r),(cid:1),φ(p),(cid:0)v∥,j,(cid:90),(cid:90),−,Φ(rp),jk,(cid:18)∂2H,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:18),v⊥,∂H,∂v⊥,(cid:19)(cid:19),v⊥,dv⊥,dv∥,=,−4π,(cid:90),(cid:90),F,v⊥,dv⊥,dv∥.,(74),Integrating,by,parts,and,neglecting,the,boundary,terms,(which,either,vanish,or,will,be,cancelled,at,the,assembly,step),we,have,that,(cid:90),(cid:90),(cid:32),∂Φ(rp),jk,∂v∥,∂H,∂v∥,+,∂Φ(rp),jk,∂v⊥,∂H,∂v⊥,(cid:33),v⊥,dv⊥,dv∥,=,−4π,(cid:90),(cid:90),F,v⊥,dv⊥,dv∥.,(75),Defining,the,matrices,and,K,(s),⊥nk,=,−,(cid:90),v(s),⊥u,v(s),⊥l,∂φ(s),n,∂v⊥,∂φ(s),k,∂v⊥,v⊥dv⊥,K,(s),∥nk,=,−,(cid:90),v(s),∥u,v(s),∥l,∂φ(s),n,∂v∥,∂φ(s),k,∂v∥,dv∥,(76),(77),Report,2070839-TN-07,and,expanding,(cid:88),(cid:88),F,=,Φ(rp),jk,F,rp,jk,,,H,=,(cid:88),(cid:88),Φ(rp),jk,H,rp,jk,,,rp,jk,rp,jk,we,find,that,the,row,of,the,unassembled,matrix,is,(cid:88),(cid:16),mn,K,(r),∥jmM,(p),⊥kn,+,K,(p),⊥knM,(r),∥jm,(cid:17),H,rp,mn,=,−4π,(cid:88),mn,M,(r),∥jmM,(p),⊥knF,rp,mn.,14,(78),(79),We,impose,Dirichlet,boundary,conditions,on,the,assembled,matrices,using,the,exact,values,of,the,required,functions.,Once,the,coefficients,H,rp,jk,are,known,then,the,same,matrices,can,be,used,in,an,identical,fashion,to,solve,for,Grp,jk.,A,similar,matrix,equation,can,be,written,down,to,solve,for,Grp,10jk,,with,the,only,difference,being,in,the,source,terms,on,the,right-hand,side.,Explicitly,,these,results,are,20jk,and,H,rp,(cid:88),(cid:16),mn,K,(r),∥jmM,(p),⊥kn,+,K,(p),⊥knM,(r),∥jm,(cid:17),H,rp,01mn,=,−4π,(cid:88),mn,∥jmM,(p),P,(r),⊥knF,rp,mn.,(80),where,we,have,defined,P,(r),∥mj,=,(cid:90),v(r),⊥u,v(r),⊥l,φ(r),m,∂φ(r),j,∂v∥,dv∥,,and,(cid:88),(cid:16),mn,K,(r),∥jmM,(p),⊥kn,+,K,(p),⊥knM,(r),∥jm,(cid:17),Grp,mn,=,2,(cid:88),mn,M,(r),∥jmM,(p),⊥knH,rp,mn.,(81),(82),Note,that,a,double,Poisson,solve,is,required,to,obtain,G,from,F,.,We,choose,to,find,∂H/∂v∥,by,a,separate,Poisson,solve,rather,than,differentiating,H,to,improve,numerical,accuracy.,To,find,the,equations,for,the,other,coefficients,in,the,fluxes,,we,must,repeat,the,exercise,above.,The,PDE,defining,G01,=,∂G/∂v⊥,is,distinct,from,Poisson’s,equation,and,will,require,different,matrix,elements,on,the,left-hand,side.,We,again,start,by,multiplying,equation,(30),by,the,function,v2,jk,,,we,find,that,⊥Φ(rp),(cid:90),(cid:90),⊥Φ(rp),v2,jk,(cid:18),∂2G01,∂v∥,2,+,1,v⊥,∂,∂v⊥,(cid:90),(cid:90),(cid:18),v⊥,∂G01,∂v⊥,(cid:19),−,G01,v2,⊥,(cid:19),dv⊥,dv∥,v2,⊥,∂H,∂v⊥,dv⊥,dv∥.,(83),=,2,Integrating,by,parts,,expanding,the,derivatives,,and,neglecting,boundary,terms,,we,find,that,(cid:90),(cid:90),(cid:32),−,v2,⊥,∂Φ(rp),jk,∂v∥,(cid:90),(cid:90),∂G01,∂v∥,∂H,∂v⊥,v2,⊥,=,2,+,v2,⊥,∂G01,∂v⊥,∂Φ(rp),jk,∂v⊥,+,v⊥Φ(rp),jk,∂G01,∂v⊥,(cid:33),+,Φ(rp),jk,G01,dv⊥,dv∥,(84),dv⊥,dv∥.,Report,2070839-TN-07,Expanding,G01,=,(cid:80),rp,(cid:80),mn,Φ(rp),mn,Grp,01mn,we,thus,find,the,unassembled,matrix,row,(cid:88),−,Grp,01mn,mn,(cid:90),(cid:90),(cid:32),v2,⊥,∂Φ(rp),jk,∂v∥,∂Φ(rp),mn,∂v∥,(cid:90),(cid:90),+,v2,⊥,∂Φ(rp),mn,∂v⊥,∂Φ(rp),jk,∂v⊥,+,v⊥Φ(rp),jk,∂Φ(rp),mn,∂v⊥,⊥Φ(rp),v2,jk,∂Φ(rp),mn,∂v⊥,dv⊥,dv∥.,+,Φ(rp),jk,Φ(rp),mn,=,2,(cid:88),mn,H,rp,mn,If,we,define,the,matrix,elements,R(s),⊥nk,=,n,(v⊥)φ(s),φ(s),k,(v⊥),v2,⊥dv⊥,,N,(s),⊥nk,=,(cid:90),v(s),⊥u,v(s),⊥l,n,(v⊥)φ(s),φ(s),k,(v⊥),dv⊥,,∂φ(s),n,∂v⊥,∂φ(s),k,∂v⊥,⊥dv⊥,,P,(s),v2,⊥nk,=,(cid:90),v(s),⊥u,v(s),⊥l,φ(s),n,∂φ(s),k,∂v⊥,v⊥dv⊥,,(cid:90),v(s),⊥u,v(s),⊥l,(cid:90),v(s),⊥u,v(s),⊥l,J,(s),⊥nk,=,−,and,U,(s),⊥nk,=,(cid:90),v(s),⊥u,v(s),⊥l,φ(s),n,∂φ(s),k,∂v⊥,v2,⊥dv⊥,then,we,can,write,the,row,as,K,(r),∥jmR(p),⊥kn,+,M,(r),∥jmJ,(p),⊥kn,−,M,(r),∥jmP,(p),⊥kn,−,M,(r),∥jmN,(p),⊥kn,(cid:17),Grp,01mn,(cid:88),(cid:16),mn,=,2,(cid:88),mn,M,(r),∥jmU,(p),⊥knH,rp,mn.,15,(cid:33),dv⊥,dv∥,(85),(86),(87),(88),The,PDE,for,G11,=,∂2G/∂v∥∂v⊥,is,identical,to,the,one,for,G01,apart,from,the,right-hand,side.,We,obtain,the,weak,form,∂G11,∂v∥,+,v2,⊥,∂G11,∂v⊥,∂Φ(rp),jk,∂v⊥,+,v⊥Φ(rp),jk,∂G11,∂v⊥,(cid:33),+,Φ(rp),jk,G11,dv⊥,dv∥,(cid:90),(cid:90),(cid:32),−,v2,⊥,∂Φ(rp),jk,∂v∥,(cid:90),(cid:90),=,2,v2,⊥,∂2H,∂v∥∂v⊥,dv⊥,dv∥.,The,corresponding,matrix,row,is,(cid:88),(cid:16),mn,K,(r),∥jmR(p),⊥kn,+,M,(r),∥jmJ,(p),⊥kn,−,M,(r),∥jmP,(p),⊥kn,−,M,(r),∥jmN,(p),⊥kn,=,2,(cid:88),mn,∥jmU,(p),P,(r),⊥knH,rp,mn.,Similarly,,the,PDE,for,H01,,equation,(39),has,the,weak,form,(cid:88),(cid:16),mn,K,(r),∥jmR(p),⊥kn,+,M,(r),∥jmJ,(p),⊥kn,−,M,(r),∥jmP,(p),⊥kn,−,M,(r),∥jmN,(p),⊥kn,(cid:17),Grp,11mn,(cid:17),H,rp,01mn,=,−4π,(cid:88),mn,M,(r),∥jmU,(p),⊥knF,rp,mn.,(89),(90),(91),Report,2070839-TN-07,16,To,complete,the,set,of,coefficients,,we,can,either,use,the,algebraic,equation,(35),2.,First,considering,equation,(35),and,or,the,elliptic,problem,(36),for,G02,=,∂2G/∂v⊥,projecting,on,to,the,basis,functions,,we,have,that,(cid:90),(cid:90),Φ(rp),jk,G02,v⊥dv⊥dv∥,=,(cid:90),(cid:90),Φ(rp),jk,(v⊥(2H,−,G20),−,G01),dv⊥dv∥.,(92),The,corresponding,matrix,row,therefore,becomes,⊥kn(2H,rp,∥jmM,(p),∥jmM,(p),⊥knGrp,02mn,=,M,(r),M,(r),(cid:88),(cid:88),mn,−,Grp,20mn),−,M,(r),∥jmN,(p),⊥knGrp,01mn.,(93),mn,mn,Secondly,,considering,equation,(36),,we,have,the,weak-form,(cid:90),(cid:90),(cid:32),−,v2,⊥,(cid:90),(cid:90),(cid:32),=,−,∂Φ(rp),jk,∂v∥,∂Φ(rp),jk,∂v⊥,v2,⊥,∂G02,∂v∥,+,v2,⊥,∂Φ(rp),jk,∂v⊥,∂G02,∂v⊥,−,v⊥Φ(rp),jk,(cid:33),−,4Φ(rp),jk,G02,dv⊥dv∥,∂G02,∂v⊥,(cid:33),∂H,∂v⊥,+,2v⊥Φ(rp),jk,∂H,∂v⊥,+,2Φ(rp),jk,H,dv⊥dv∥,+,(cid:90),(cid:90),Φ(rp),jk,G20,dv⊥dv∥,,(94),where,we,have,carried,out,the,integration,by,parts,and,neglected,to,write,boundary,terms.,The,corresponding,row,of,the,unassembled,matrix,therefore,becomes,⊥kn,−,P,(p),⊥kn,−,4N,(p),⊥kn,+,M,(r),(cid:17),⊥kn),∥jm(J,(p),∥jmR(p),K,(r),Grp,(cid:88),02mn,(cid:16),(cid:32),mn,=,2,(cid:88),mn,M,(r),∥jm(J,(p),⊥kn,−,2P,(p),⊥kn,−,N,(p),⊥kn)H,rp,mn,−,2,(cid:88),mn,M,(r),∥jmN,(p),⊥knGrp,20mn,(95),(cid:33),.,3.6.,Velocity,space,integration,in,the,spectral,element,scheme,To,compute,the,boundary,data,for,the,elliptic,solvers,obtained,in,the,last,section,,we,need,to,integrate,a,function,F,=,F,(v′,⊥),multiplied,by,a,kernel,function,∥,,v′,G,=,G(v∥,,v⊥,,v′,∥,,v′,⊥),,i.e.,,we,wish,to,compute,I,=,(cid:90),∞,(cid:90),∞,−∞,0,G(v∥,,v⊥,,v′,∥,,v′,⊥)F,(v′,∥,,v′,⊥)v′,⊥dv′,⊥dv′,∥.,(96),(97),We,expand,F,in,the,Lagrange,polynomial,basis,functions,using,equation,(62),and,thus,obtain,that,(cid:88),(cid:88),I,=,jk,I,(rp),F,rp,jk,rp,with,the,integration,over,local,elements,jk,I,(rp),jk,=,(cid:90),v(r),∥u,(cid:90),v(p),⊥u,v(r),∥l,v(p),⊥l,G(v∥,,v⊥,,v′,∥,,v′,⊥)φ(r),j,(v′,∥)φ(p),k,(v′,⊥)v′,⊥dv′,⊥dv′,∥.,(98),(99),When,assembling,the,element,problem,,one,must,recall,that,at,element,boundaries,the,nodal,value,of,F,has,an,interpolating,polynomial,that,should,be,integrated,over,on,both,elements.,Report,2070839-TN-07,17,4.,Boundary,conditions,for,time,evolution,The,collision,operator,described,in,this,report,is,intended,to,be,used,alongside,explicit,pseudo-spectral,advection,terms,in,both,space,and,velocities.,As,such,,boundary,conditions,are,imposed,explicitly,at,the,boundaries,of,the,domain,each,time,step,,by,explicitly,forcing,the,appropriate,conditions,on,Fs.,In,the,assembly,of,the,collision,operator,(71),,we,imposed,a,“zero-flux”,boundary,condition,Γ∥(v∥,=,±L∥),=,Γ⊥(v⊥,=,L⊥),=,0,by,virtue,of,neglecting,the,boundary,fluxes,in,the,implementation.,However,,for,the,advection,terms,,we,must,impose,Fs,=,0,at,the,upwind,boundaries,of,velocity,space.,For,physical,solutions,we,would,expect,Fs,=,0,and,the,“zero-flux”,boundary,conditions,to,be,equivalent.,The,use,of,cylindrical,coordinates,introduces,a,spurious,boundary,into,the,coordinate,system,at,v⊥,=,0.,To,ensure,regularity,of,the,solution,at,v⊥,=,0,,we,impose,(cid:12),(cid:12),(cid:12),(cid:12)v⊥=0,We,motivate,this,by,noting,that,for,Fs,=,Fs(v),to,be,continuous,and,differentiable,as,a,function,of,the,vector,velocity,argument,v,near,v,=,0,,then,(100),must,be,satisfied.,∂Fs,∂v⊥,(100),=,0.,5.,Ad-hoc,numerical,conserving,terms,The,numerical,scheme,here,is,chosen,for,speed,and,accuracy,rather,than,for,exact,satisfaction,of,the,necessary,conservative,properties,(5)-(7).,To,ensure,that,the,numerical,scheme,also,preserves,the,density,,parallel,velocity,and,pressure,at,each,for,time-evolving,simulations,we,introduce,ad-,time,step,(to,machine,precision),,hoc,conserving,terms,which,make,a,correction,which,is,(in,principle),of,order,the,discretisation,error.,Noting,the,definitions,of,the,plasma,density,ns,and,parallel,flow,u∥,s,,(cid:90),(cid:90),ns,=,Fs,2πv⊥dv⊥,dv∥,,and,respectively,,we,define,(cid:90),(cid:90),nsu∥,s,=,v∥Fs,2πv⊥dv⊥,dv∥,,Css,[Fs,,Fs],=,C,∗,ss,[Fs,,Fs],−,(cid:0)x0,+,x1(v∥,−,u∥,s),+,x2,(cid:0)(v∥,−,u∥,s)2,+,v2,⊥,(cid:1)(cid:1))Fs,,(103),where,C,∗,ss,[Fs,,Fs],denotes,the,numerically,calculated,finite-element,collision,operator,given,by,∂F/∂t,in,equation,(71),,and,the,coefficients,x0,,x1,and,x2,are,determined,by,the,requirements,that,(5)-(7),are,exactly,satisfied.,These,requirements,lead,to,the,matrix,equation,,,,,,,,,ns,0,3ps,0,p∥,s,q∥,s,3ps,q∥,s,˜Rs,x0,x1,x2,,,,,,,=,,,∆ns,ns∆u∥,s,−,u∥,s∆ns,3∆ps,,,(104),(101),(102),(106),(107),Report,2070839-TN-07,18,where,the,vector,components,on,the,right,hand,side,are,the,moments,of,C,∗,ss,[Fs,,Fs],(cid:90),(cid:90),∆ns,=,C,∗,ss,[Fs,,Fs],2πv⊥dv⊥dv∥,,∆u∥,s,=,∆ps,=,1,ns,1,3,(cid:90),(cid:90),(cid:90),(cid:90),v∥C,∗,ss,[Fs,,Fs],2πv⊥dv⊥dv∥,,and,(105),(cid:0)(v∥,−,u∥,s)2,+,v2,⊥,(cid:1),C,∗,ss,[Fs,,Fs],2πv⊥dv⊥dv∥,,and,the,compoments,of,the,matrix,on,the,left-hand,side,are,given,by,the,moments,of,Fs.,We,have,that,the,total,pressure,ps,=,(2p⊥,s,+,p∥,s)/3,with,the,parallel,and,perpendicular,pressures,given,by,(cid:90),(cid:90),p∥,s,=,(v∥,−,u∥,s)2Fs,2πv⊥dv⊥,dv∥,,and,p⊥,s,=,(cid:90),(cid:90),1,2,v2,⊥Fs,2πv⊥dv⊥,dv∥,,respectively.,The,parallel,heat,flux,is,given,by,(cid:90),(cid:90),q∥,s,=,(v∥,−,u∥,s)((v∥,−,u∥,s)2,+,v2,⊥)Fs,2πv⊥dv⊥,dv∥,,(108),and,the,higher-order,moment,(cid:90),(cid:90),˜Rs,=,((v∥,−,u∥,s)2,+,v2,⊥)2Fs,2πv⊥dv⊥,dv∥.,(109),To,achieve,good,results,with,this,method,,the,boundary,conditions,that,are,imposed,on,Fs,should,also,be,imposed,on,C,∗,ss,[Fs,,Fs],before,evaluating,the,values,of,the,set,of,coefficients,{xi}.,We,note,the,similarity,of,these,ad-hoc,conserving,terms,to,those,employed,for,similar,reasons,where,the,density,,parallel,flow,,and,pressure,are,required,to,be,conserved,exactly,[10,,21].,6.,Numerical,implementation,and,results,We,have,implemented,an,explicit,form,of,the,weak-form,collision,operator,in,a,script,to,test,the,numerical,properties,,as,well,as,in,the,main,time-advance,loop,of,“moment,kinetics”.,Specifically,,we,have,implemented,the,assembled,weak-,form,problems,defined,by,equations,(71),,(79),,(80),,(82),,(88),,(90),,(91),,(95),,using,sparse,matrices,https://docs.julialang.org/en/v1/stdlib/SparseArrays/,,using,the,“moment,kinetics”,shared-memory,MPI,,with,appropriate,calculations,of,the,boundary,data,using,the,integration,weights,defined,by,(99).,We,use,Gauss-Legendre,polynomials,to,define,the,Lobatto,and,Radau,collocation,grid,points.,We,have,implemented,the,scheme,for,arbitrary,order,of,polynomial.,The,implementation,may,be,found,in,the,branch,at,the,following,URL,https://,Report,2070839-TN-07,19,github.com/mabarnes/moment_kinetics/tree/merge_fkpl_collisions.,The,soft-,ware,required,to,execute,the,tests,described,in,section,6.1,are,contained,in,the,test,script,in,the,following,URL,https://github.com/mabarnes/moment_kinetics/,tree/merge_fkpl_collisions/test_scripts/2D_FEM_assembly_test.jl,,on,com-,mit,c24dbc6.,relaxation,test,described,in,section,6.2,uses,an,input,file,based,on,https://github.com/mabarnes/moment_kinetics/blob/merge_fkpl_,collisions/examples/fokker-planck/fokker-planck-relaxation.toml.,This,6.1.,Evaluation,tests,We,wish,to,test,the,three,properties,of,the,collision,operator,(5)-(7).,To,facilitate,this,test,we,define,three,quantities,which,measure,the,change,in,the,moments,of,the,distribution,function,due,to,the,collision,operator,,given,by,equation,(105).,We,test,whether,or,not,the,collision,operator,vanishes,on,a,prescribed,Maxwellian,distribution,,i.e.,,Css,(cid:2)F,M,s,,,F,M,s,(cid:3),=,0,,(110),as,a,proxy,for,testing,the,entropy,production,inequality,(8).,In,figures,1,,2,,3,,and,4,,we,carry,out,the,test,on,2,cores,for,varying,Nelement,at,fixed,Ngrid,=,3,,5,,7,,and,9,,respectively.,Here,Nelement,is,the,number,of,elements,in,the,v⊥,dimension,and,half,the,number,of,elements,in,the,v∥,dimension.,The,quantity,Ngrid,is,the,number,of,points,per,element.,We,take,the,maximum,velocity,to,be,L∥,=,L⊥,=,6.0.,Note,that,reducing,the,maximum,v∥,and,v⊥,on,the,grids,for,a,fixed,integrand,reduces,the,accuracy,of,the,numerical,integration,because,the,true,velocity,integrals,should,extend,to,infinite,velocities.,We,choose,to,carry,out,the,test,with,a,Maxwellian,with,a,normalised,density,ns/nref,=,1.0,,a,normalised,u∥,s/cref,=,1.0,,with,cref,=,(cid:112)2Tref/mref,and,a,normalised,Ts/Tref,=,1.0.,In,the,figure,(a),of,figures,1-4,,we,plot,both,the,infinity,norm,of,the,error,ϵ∞,and,the,L2,norm,of,the,error,ϵL2,of,calculating,the,collision,operator,with,respect,to,the,expected,value,(which,is,zero).,We,see,that,the,infinity,norm,gives,a,larger,value,than,the,L2,norm,in,all,cases,by,a,factor,of,an,order,of,magnitude.,This,is,due,to,numerical,oscillations,near,v⊥,=,0,where,the,differential,equations,become,singular.,For,Ngrid,=,3,(quadratic,interpolating,polynomials),the,error,in,computing,the,collision,operator,does,not,decrease,rapidly,with,Nelement,,but,the,error,does,follow,the,expected,scaling,for,differentiation,(cid:19)Ngrid−1,(cid:18),1,Nelement,.,(111),However,,the,quantities,∆ns,,∆u∥,s,,and,∆ps,approach,zero,rapidly,at,(or,better,than),the,expected,scaling,for,numerical,integration,errors,(cid:18),1,Nelement,(cid:19)Ngrid+1,.,(112),This,picture,becomes,clearer,for,larger,Ngrid,,as,evidenced,in,figures,2-4.,Report,2070839-TN-07,20,To,demonstrate,the,attained,performance,of,the,explicit,collision,operator,,in,figure,(b),of,figures,1-4,,we,plot,the,timing,data,(in,milli,seconds),for,completing,the,initialisation,and,evaluation,of,the,collision,operator.,The,expected,scaling,for,the,initialisation,is,N,3,element,,by,virtue,of,the,calculation,of,the,integration,weights,for,the,boundary,data.,The,expected,scaling,for,the,evaluation,of,the,collision,operator,depends,on,which,operation,dominates,the,calculation.,If,it,is,the,computation,of,the,boundary,data,it,is,N,3,element,,whereas,if,it,is,the,elliptic,solve,or,the,assembly,of,the,right,hand,side,of,equation,(71),then,the,scaling,would,be,expected,to,be,N,2,element,due,to,the,sparse,nature,of,these,operations.,In,figure,1b,for,Ngrid,=,3,the,timing,for,both,the,initialisation,and,single,evaluation,scales,like,N,3,element,,whereas,,in,figures,2b-4b,for,Ngrid,=,5-9,,respectively,,we,see,that,a,scaling,closer,to,N,2,element,is,acheived,for,the,evaluation,step.,This,is,promising,for,the,scaling,of,the,operator,to,large,problem,sizes.,To,understand,the,dominant,source,of,the,numerical,error,,we,find,it,useful,to,plot,the,infinity,and,L2,norm,error,measures,of,the,numerically,calculated,2,,∂G/∂v⊥,,∂2G/∂v∥∂v⊥,,Rosenbluth,potential,coefficients,∂H/∂v∥,,∂H/∂v⊥,,∂2G/∂v∥,and,∂2G/∂v⊥,2.,The,exact,values,are,known,and,may,be,computed,easily,for,shifted,Maxwellian,distributions,[22].,We,plot,this,data,for,Ngrid,=,3,in,figure,5,,and,for,Ngrid,=,9,in,in,figure,6.,We,see,that,for,both,Ngrid,=,3,and,Ngrid,=,9,,the,L2,norm,error,is,smaller,by,one,or,two,orders,of,magnitude,than,the,infinity,norm,error.,This,is,due,to,numerical,oscillations,near,v⊥,=,0.,However,,in,both,cases,the,errors,decay,to,zero,at,least,at,the,rate,given,by,equation,(111),(for,Ngrid,=,3),,or,at,the,proper,expected,rate,given,by,equation,(112),(for,Ngrid,=,9).,Note,that,our,numerical,calculation,of,the,boundary,data,does,involve,a,numerical,differentiation,of,F,,,see,equations,(47)-(52).,In,the,case,of,Ngrid,=,9,,we,see,that,the,errors,deviate,from,the,scaling,for,high,resolution,,which,we,take,to,mean,that,we,have,reached,the,smallest,error,possible,with,double,floating,point,precision.,We,take,this,to,indicate,that,the,dominant,source,of,numerical,error,in,evaluating,the,collision,operator,in,fact,comes,from,the,unavoidable,numerical,differentiation,,rather,than,from,the,errors,in,obtaining,the,Rosenbluth,potential,coefficients.,Indeed,,comparable,levels,of,error,to,that,seen,in,computing,the,collision,operator,may,be,obtained,by,simply,using,the,weak,method,to,differentiate,F,to,find,the,second,derivatives,in,v∥,and,v⊥.,6.2.,Relaxation,to,a,Maxwellian,distribution:,testing,equation,(8),It,is,important,to,test,whether,or,not,the,numerical,self-collision,operator,can,provide,a,stable,,steady-state,numerical,solution,which,is,close,to,a,Maxwellian,distribution,,whilst,satisfying,equation,(8).,Using,the,ad-hoc,numerically,conserving,model,(103),,we,find,that,we,can,obtain,a,stable,solution,even,for,very,low,numerical,resolution.,In,figure,7,,we,show,time,traces,of,the,change,in,the,density,,parallel,flow,and,pressure,over,the,course,of,the,simulation.,In,figure,8,we,show,the,entropy,production,,calculated,using,the,definition,(8),and,using,the,following,approximation,for,the,logarithm,of,the,Report,2070839-TN-07,21,(a),(b),Figure,1:,The,numerical,error,and,timing,data,for,the,test,carried,out,on,2,cores,with,Ngrid,=,ng,=,3,points,per,element.,The,infinity,and,L2,norms,of,the,collision,operator,are,shown,and,compared,to,the,expected,scalings,for,differentiation,and,integration,(111),and,(112),,respectively.,The,timing,data,for,the,initialisation,(init),and,a,single,evaluation,of,the,collision,operator,(step),is,given,in,milliseconds.,Report,2070839-TN-07,22,(a),(b),Figure,2:,The,numerical,error,and,timing,data,for,the,test,carried,out,on,2,cores,with,Ngrid,=,ng,=,5,points,per,element.,Report,2070839-TN-07,23,(a),(b),Figure,3:,The,numerical,error,and,timing,data,for,the,test,carried,out,on,2,cores,with,Ngrid,=,ng,=,7,points,per,element.,Report,2070839-TN-07,24,(a),(b),Figure,4:,The,numerical,error,and,timing,data,for,the,test,carried,out,on,2,cores,with,Ngrid,=,ng,=,9,points,per,element.,Report,2070839-TN-07,25,(a),(b),Figure,5:,In,figure,(a),we,plot,the,infinity,norm,of,the,error,ϵ∞,in,computing,the,Rosenbluth,potential,coefficients,,for,Ngrid,=,3.,In,figure,(b),we,plot,the,L2,norm,of,the,error,ϵL2,for,Ngrid,=,3.,Report,2070839-TN-07,26,(a),(b),Figure,6:,In,figure,(a),we,plot,the,infinity,norm,of,the,error,ϵ∞,in,computing,the,Rosenbluth,potential,coefficients,,for,Ngrid,=,9.,In,figure,(b),we,plot,the,L2,norm,of,the,error,ϵL2,for,Ngrid,=,9,Report,2070839-TN-07,27,Figure,7:,The,changes,in,the,first,three,moments,of,the,distribution,function,ns,,u∥,s,,and,ps,as,a,result,of,time,evolution,with,the,Fokker-Planck,collision,operator,defined,by,equations,(71),and,(103),(i.e.,,with,the,numerical,conserving,terms).,The,moments,are,well,conserved,,despite,the,low,resolution,used:,here,we,use,Ngrid,=,5,,Nelement,=,4,elements,in,the,v⊥,dimension,and,2Nelement,=,8,elements,in,the,v∥,dimension.,We,take,the,maximum,velocity,to,be,L∥,=,L⊥,=,6.0.,We,take,∆t,=,10−3.,distribution,function,ln,F,=,(cid:88),(cid:88),rp,ij,ln,(cid:0)|F,rp,jk,|,+,ϵ(cid:1),Φ(rp),jk,,,(113),sc3,where,ϵ,=,10−15.,The,approximation,(113),is,adequate,if,the,solution,is,converging,with,increasing,resolution,in,a,strong,sense.,In,figure,9,we,show,the,L2,norm,of,Fs,−,F,M,s,.,The,simulation,uses,a,collision,frequency,νss,=,γssnref/m2,ref,=,cref/Lref,,and,is,run,for,a,time,of,200Lref/cref.,We,use,Ngrid,=,5,,Nelement,=,4,elements,in,the,v⊥,domain,,2Nelement,=,8,elements,in,the,v∥,domain,,and,we,take,the,maximum,velocity,to,be,L∥,=,L⊥,=,3.0.,We,take,∆t,=,10−3,and,we,use,the,RK4,explicit,time,integration,method.,Despite,the,low,resolution,,the,numerical,solution,is,stable.,The,small,errors,in,the,moments,are,of,order,10−7,at,t,=,200Lref/cref,,which,is,larger,than,machine,precision.,To,reach,this,time,,200,×,1000,×,4,≃,108,evaluations,of,Css,[Fs,,Fs],have,been,carried,out.,This,might,explain,the,errors,in,the,moments,if,the,correction,terms,have,a,machine-precision,≃,10−15,systematic,bias.,Simulations,with,increasing,numerical,resolution,show,a,steady,states,with,smaller,L2(Fs,−,F,M,s,).,Report,2070839-TN-07,28,Figure,8:,We,plot,the,the,entropy,production,˙S,,defined,in,equation,(8).,Note,that,˙S,remains,positive,and,tends,to,0+.,The,resolutions,are,provided,in,the,main,text.,Figure,9:,We,plot,the,L2,norm,of,Fs,−,F,M,increasingly,close,to,F,M,The,resolutions,are,provided,in,the,main,text.,s,.,This,figure,indicates,that,Fs,becomes,s,before,converging,on,a,steady-state,numerical,Maxwellian.,Report,2070839-TN-07,7.,Discussion,and,outlook,29,In,this,report,we,have,investigated,a,particular,weak-form,representation,of,the,explicit,Landau,Fokker-Planck,collision,operator.,We,choose,to,use,the,RMJ,form,of,the,Fokker-,Planck,operator,to,permit,the,use,of,sparse,elliptic,solves,for,determining,the,coefficients,of,the,nonlinear,operator.,We,have,demonstrated,that,this,choice,can,lead,to,an,optimal,scaling,of,the,cost,of,evaluating,the,operator,for,a,single,time,step,∝,N,2,element,with,the,number,of,elements,Nelement,in,each,of,the,velocity,space,dimensions,v∥,and,v⊥.,We,also,demonstrated,a,successful,time-evolving,simulation,with,low,resolution,,demonstrating,that,the,self-collision,operator,can,successfully,relax,the,distribution,function,to,a,stable,steady,state,that,is,close,to,a,Maxwellian,distribution.,Unfortunately,,the,order,unity,factor,in,front,of,the,expected,cost,scaling,is,large,for,all,Ngrid.,To,permit,the,collision,operator,to,be,routinely,used,in,the,time,evolving,“moment,kinetics”,code,alongside,other,physics,features,we,would,ideally,want,to,further,optimise,the,implementation,for,speed.,This,might,be,achieved,with,an,improvement,to,the,shared-memory,paralellism,,distributed-memory,parallelisation,across,nodes,or,through,optimisations,of,the,numerical,method,for,speed,at,the,cost,of,accuracy.,The,use,of,distributed,memory,to,parallelise,the,collision,operator,calculation,should,be,fairly,straightforward.,The,dominant,costs,which,contribute,to,the,time,taken,to,evaluate,the,operator,are,the,calculation,of,the,boundary,data,and,the,assembly,of,the,right-hand,side,of,equation,(71).,Both,of,these,steps,are,embarrassingly,parallel.,Possible,optimisations,of,the,numerical,method,could,involve,the,following.,First,,choosing,to,determine,the,boundary,data,for,the,elliptic,solves,using,a,multipole,expansion,of,the,Green’s,functions,definition,of,G,and,H,,equations,(16),and,(17),,respectively.,This,method,may,permit,the,evaluation,of,the,boundary,data,using,only,an,order,unity,number,of,velocity,integrals,,providing,the,maximum,value,of,v∥,and,v⊥,on,the,grid,,L∥,and,L⊥,,respectively,,are,sufficiently,large.,Second,,choosing,to,evaluate,the,boundary,data,at,fewer,locations,and,constructing,a,larger-scale,interpolation,of,the,coefficients,on,the,boundary,might,save,computation,time,without,sacrificing,significant,accuracy,,again,if,L∥,and,L⊥,are,large,enough,for,the,Rosenbluth,potential,coefficients,to,have,a,simple,functional,form.,Finally,,one,could,choose,to,use,a,different,interpolation,scheme,defining,the,right-hand,side,of,equation,(71).,One,could,consider,using,the,commonly,used,spectral-element,method,“quadrature,crime”,of,assuming,that,mass,matrices,are,diagonal,[23],to,reduce,the,number,of,operations,due,to,the,nonlinear,stiffness,matrices,defined,by,equations,(72),and,(73).,Appendix,A.,Evaluating,the,gyroaveraged,functions,To,see,how,to,evaluate,the,required,gyroaveraged,functions,IG1,,IH0,,IH1,,and,IH2,,consider,IG0,=,1,2π,(cid:90),π,−π,g,dϑ′,,(A.1),Report,2070839-TN-07,as,well,as,and,IG2,=,1,2π,(cid:90),π,−π,g,(e⊥,·,e′,⊥)2,dϑ′,,IG3,=,1,2π,(cid:90),π,−π,g,(e⊥,·,e′,⊥,×,b)2,dϑ′.,30,(A.2),(A.3),We,note,that,e⊥,·,e′,that,⊥,=,cos(ϑ′,−,ϑ),and,e⊥,·,e′,⊥,×,b,=,sin(ϑ,−,ϑ′).,Expanding,g,,we,have,IG0(v∥,,v⊥,,v′,∥,,v′,⊥),=,1,2π,(cid:90),π,−π,(cid:16)(cid:0)v∥,−,v′,∥,(cid:1)2,+,v2,⊥,+,v′,⊥,2,−,2v⊥v′,⊥,cos,(ϑ′,−,ϑ),(cid:17)1/2,dϑ′.,(A.4),Here,we,can,recognise,an,Elliptic,integral.,Suitable,rearrangement,and,relabeling,give,us,IG0(v∥,,v⊥,,v′,∥,,v′,⊥),=,(cid:16)(cid:0)v∥,−,v′,∥,(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)1/2,2,π,E(m(v∥,,v⊥,,v′,∥,,v′,⊥)),with,m(v∥,,v⊥,,v′,∥,,v′,⊥),=,4v⊥v′,⊥,(cid:16)(cid:0)v∥,−,v′,∥,(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)−1,and,we,have,used,the,definition,of,the,complete,elliptic,integral,of,the,first,kind,K(m),=,(cid:90),π/2,0,1,1,−,m,sin2,θ,(cid:112),dθ,Complete,elliptic,integral,of,the,second,kind,E(m),=,(cid:90),π/2,0,(cid:112),1,−,m,sin2,θ,dθ,The,remaining,integrals,are,(A.5),(A.6),(A.7),(A.8),IG1(v∥,,v⊥,,v′,∥,,v′,⊥),=,−,2,π,(cid:16)(cid:0)v∥,−,v′,∥,(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)1/2,(cid:18),2,−,m,3m,E(m),−,2,3,(1,−,m),m,(cid:19),K(m),(A.9),IG2(v∥,,v⊥,,v′,∥,,v′,IG3(v∥,,v⊥,,v′,∥,,v′,⊥),=,2,π,⊥),=,2,π,(cid:16)(cid:0)v∥,−,v′,∥,(cid:16)(cid:0)v∥,−,v′,∥,(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)1/2,1,15m2,(cid:0)(7m2,+,8m,−,8)E(m),+,4(2,−,m)(1,−,m)K(m)(cid:1),(A.10),(cid:1)2,+,(v⊥,+,v′,⊥)2(cid:17)1/2,1,15m2,(cid:0)8(m2,−,m,+,1)E(m),−,4(2,−,m)(1,−,m)K(m)(cid:1),(A.11),Report,2070839-TN-07,31,where,we,have,used,the,identities,(cid:90),π/2,0,sin2,θ,(cid:112),1,−,m,sin2,θ,dθ,=,1,−,m,3m,K(m),+,2m,−,1,3m,E(m),(A.12),(cid:90),π/2,(cid:112),sin4,θ,0,1,−,m,sin2,θ,dθ,=,1,15m2,(cid:0)2(2m,+,1)(1,−,m)K(m),+,(8m2,−,3m,−,2)E(m)(cid:1),(cid:90),π/2,0,(1,−,2,sin2,θ),(cid:112),1,−,m,sin2,θ,dθ,=,2,−,m,3m,E(m),−,2,3,(1,−,m),m,K(m),(cid:0)(7m2,+,8m,−,8)E(m),+,4(2,−,m)(1,−,m)K(m)(cid:1),(A.13),(A.14),(A.15),(cid:90),π/2,0,(1−2,sin2,θ)2(cid:112),1,−,m,sin2,θ,dθ,=,1,15m2,(cid:90),π/2,0,(cid:0)1,−,(1,−,2,sin2,θ)2(cid:1)(cid:112),1,−,m,sin2,θ,dθ,=,1,15m2,(cid:0)8(m2,−,m,+,1)E(m),−,4(2,−,m)(1,−,m)K(m)(cid:1),(A.16),If,we,wish,to,evaluate,directly,the,integral,equations,for,∂Hs′/∂v∥,and,∂Hs′/∂v⊥,then,we,require,the,following,elliptic,integrals,(cid:90),π,IH0,=,IH1,=,1,2π,1,2π,(cid:90),π,−π,1,g,dϑ′,,−π,e⊥,·,e′,⊥,g,dϑ′,,and,IH2,=,1,2π,(cid:90),π,−π,⊥)2,(e⊥,·,e′,g,dϑ′.,Using,the,methods,described,above,,we,find,that,IH0,=,2,π,(cid:0)(v∥,−,v′,∥)2,+,(v⊥,+,v′,⊥)2(cid:1)−1/2,K(m),,(A.17),(A.18),(A.19),(A.20),IH1,=,−,2,π,(cid:0)(v∥,−,v′,∥)2,+,(v⊥,+,v′,⊥)2(cid:1)−1/2,(cid:18),m,−,2,m,K(m),+,2,m,(cid:19),E(m),,,(A.21),and,IH2,=,2,π,(cid:0)(v∥,−,v′,∥)2,+,(v⊥,+,v′,⊥)2(cid:1)−1/2,(cid:18),3m2,−,8m,+,8,3m2,K(m),+,(cid:19),4m,−,8,3m2,E(m),.,(A.22),Here,we,have,used,that,(cid:90),π/2,0,(cid:0)1,−,2,sin2,θ(cid:1),(cid:0)1,−,m,sin2,θ(cid:1)−1/2,dθ,=,m,−,2,m,K(m),+,2,m,E(m),,(A.23),and,(cid:90),π/2,0,(cid:0)1,−,2,sin2,θ(cid:1)2,(cid:0)1,−,m,sin2,θ(cid:1)−1/2,dθ,=,3m2,−,8m,+,8,3m2,K(m),+,4m,−,8,3m2,E(m).,(A.24),Report,2070839-TN-07,32,Figure,B1:,We,plot,the,infinity,norm,of,the,errors,ϵ,of,the,potentials,∂H/∂v∥,,∂H/∂v⊥,,∂G/∂v⊥,,∂2G/∂v⊥,2,,∂2G/∂v⊥∂v∥,for,a,Maxwellian,input,distribution,,compared,to,the,expected,scalings,for,differentiation,and,integration,,equations,(111),and,(112),,respectively.,2,,∂2G/∂v∥,Appendix,B.,Computing,the,Rosenbluth,potentials,by,direct,integration,A,more,direct,,but,less,efficient,,method,for,computing,the,Rosenbluth,potentials,is,to,use,the,integral,expressions,(47)-(52),for,all,(v∥,,v⊥),rather,than,just,the,boundary,values.,Here,we,show,the,results,of,such,a,calculation,to,demonstrate,the,correct,implementation,of,(47)-(52),and,the,results,in,Appendix,A.,In,figure,B1,we,plot,the,infinity-norm,errors,on,the,calculation,by,direct,integration,of,the,derivatives,of,the,Rosebluth,potentials,for,a,Maxwellian,input,,for,which,the,results,are,known,analytically,(see,e.g.,[22]).,We,see,that,the,integration,error,becomes,small,for,increasing,resolution,,indicating,that,the,definitions,of,the,integrands,are,correct.,However,,the,errors,eventually,deviate,from,the,expected,scaling.,This,is,due,to,problems,carrying,out,the,integral,accurately,in,the,region,on,the,integrand,where,v′,is,such,that,Fs(v′),∼,O,(1).,This,problem,might,be,addressed,with,an,improved,integration,quadrature,,or,by,using,higher,than,double,precision,to,compute,the,integrand.,Note,that,this,difficulty,does,not,affect,the,integration,of,the,potentials,in,the,far-field,region,at,the,velocity,space,boundary,–,meaning,that,machine-precision,accuracy,can,be,achieved,in,the,numerical,method,presented,in,the,main,text.,This,is,evident,from,figures,5,and,6.,The,script,used,to,generate,figure,B1,may,be,found,at,Report,2070839-TN-07,33,the,following,URL:,https://github.com/mabarnes/moment_kinetics/blob/merge_,fkpl_collisions/test_scripts/fkpl_direct_integration_test.jl.,[1],Kardar,M,2007,Statistical,Physics,of,Particles,(Cambridge,University,Press),[2],Dellar,P,Kinetic,Theory,University,of,Oxford,Master,Course,in,Mathematical,and,Theoretical,Physics:,lecture,series.,URL,https://people.maths.ox.ac.uk/dellar/MMPkinetic.html,[3],Rosenbluth,M,N,,MacDonald,W,M,and,Judd,D,L,1957,Phys.,Rev.,107,1–6,[4],Hazeltine,R,D,and,Meiss,J,D,2003,Plasma,Confinement,(New,York:,Dover),[5],Helander,P,and,Sigmar,D,J,2002,Collisional,Transport,in,Magnetized,Plasmas,(Cambridge,,UK:,Cambrige,University,Press),[6],Parra,F,I,Collisional,Plasma,Physics:,(I),Fokker-Planck,collision,operator,University,of,Oxford,Master,Course,in,Mathematical,and,Theoretical,Physics:,2019,lecture,series.,URL,https://www-thphys.physics.ox.ac.uk/people/FelixParra/CollisionalPlasmaPhysics/,CollisionalPlasmaPhysics.html,[7],Alouani-Bibi,F,,Shoucri,M,and,Matte,J,P,2004,Computer,physics,communications,164,60–66,[8],Pataki,A,and,Greengard,L,2011,Journal,of,Computational,Physics,230,7840–7852,[9],Hirvijoki,E,and,Adams,M,F,2017,Phys.,Plasmas,24,032121,[10],Abazorius,M,2023,University,of,Oxford,DPhil,Thesis,(In,Progress),[11],Wilkie,G,J,,Keßler,T,and,Rjasanow,S,2023,Comput.,Phys.,Commun.,291,108812,[12],Catto,P,J,and,Simakov,A,N,2004,Phys.,Plasmas,11,90,[13],Hinton,F,L,and,Hazeltine,R,D,1976,Rev.,Mod.,Phys.,48,239–308,[14],Catto,P,J,and,Tsang,K,T,1977,Phys.,Fluids,20,396–401,[15],Abel,I,G,,Barnes,M,,Cowley,S,C,,Dorland,W,and,Schekochihin,A,A,2008,Phys.,Plasmas,15,122509,[16],Barnes,M,,Abel,I,,Dorland,W,,Ernst,D,,Hammett,G,,Ricci,P,,Rogers,B,,Schekochihin,A,and,Tatsuno,T,2009,Phys.,Plasmas,16,[17],Sugama,H,,Watanabe,T,H,and,Nunami,M,2009,Phys.,Plasmas,16,112503,[18],Abel,I,G,,Plunk,G,G,,Wang,E,,Barnes,M,,Cowley,S,C,,Dorland,W,and,Schekochihin,A,A,Rep.,Prog.,Phys.,76,116201,[19],Parra,F,I,Collisionless,Plasma,Physics:,Course,in,Mathematical,and,Theoretical,Physics:,//www-thphys.physics.ox.ac.uk/people/FelixParra/CollisionlessPlasmaPhysics/,CollisionlessPlasmaPhysics.html,(II),Drift-Kinetics,University,of,Oxford,Master,2019,lecture,series.,URL,https:,[20],Boyd,J,P,2001,Chebyshev,and,Fourier,Spectral,Methods,(Dover),[21],Barnes,M,,Parra,F,I,,Hardman,M,R,and,Omotani,J,2021,Excalibur/Neptune,Report,2047357–,TN–D2.2+M2.5,[22],Hardman,M,R,,Omotani,J,,Barnes,M,,Newton,S,L,and,Parra,F,I,2023,Excalibur/Neptune,Report,2070839–TN–06,[23],Teukolsky,S,A,2015,J.,Comput.,Phys.,283,408–413

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