TN-02_Task11EllipticSolverTests
===============================

.. meta::
   :description: technical note

   :keywords: Excalibur-Neptune,report,2047356-TN-02-2,Task,1.1,Elliptic,solver,tests,Ben,Dudson,,Peter,Hill,,Ed,Higgins,,David,Dickinson,,and,Steven,Wright,University,of,York,David,Moxey,University,of,Exeter,June,16,,2021,Contents,1,Executive,summary,2,Elliptic,solvers,in,plasma,equations,2.1,Embedded,domain,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,3,Tests,of,elliptic,solvers,3.1,Correctness,tests,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,3.2,Performance,tests,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,4,Tests,in,slab,geometries,4.1,Slab,with,a,pole,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,5,Tests,in,tokamak,geometries,2,2,4,5,5,6,7,8,8,5.1,Tokamak,geometries,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,10,5.1.1,Numerical,tokamak,equilibria,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,10,5.1.2,Analytic,tokamak,equilibria,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,10,6,Time-evolving,systems,11,6.1,Shear,Alfv`en,wave,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,11,6.1.1,Energy,conservation,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,12,6.2,Geodesic,Acoustic,Wave,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,14,6.2.1,Simplified,model,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,14,6.2.2,Analytic,solution,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,16,6.2.3,Mass,and,energy,conservation,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,.,17,7,Conclusions,19,1,8,References,20,1,Executive,summary,This,report,gives,a,brief,overview,of,the,origin,of,elliptic,problems,in,plasma,physics,models,,and,the,tokamak,geometry,they,are,solved,in.,Drawing,mainly,on,experience,with,BOUT++,,approaches,to,testing,of,both,correctness,and,performance,,and,the,pros,and,cons,of,them,are,discussed.,A,series,of,geometries,which,can,be,used,for,testing,are,described,,starting,with,slabs,of,increasing,complexity,,and,then,tokamak,geometries.,Finally,two,simple,time-dependent,sets,of,equations,are,suggested,,one,for,the,shear,Alfv`en,wave,,and,the,other,the,Geodesic,Acoustic,Wave,(GAM).,These,provide,ways,to,test,the,long-term,numerical,stability,of,the,elliptic,solver,,together,with,the,boundary,conditions,parallel,and,perpendicular,to,the,magnetic,field.,2,Elliptic,solvers,in,plasma,equations,Elliptic,problems,appear,in,plasma,equations,typically,in,the,electromagnetic,fields,,in,the,limit,where,the,displacement,current,is,neglected,and,light,speed,goes,to,infinity.,Two,common,components,are,solving,for,the,electrostatic,potential,φ,,and,the,parallel,component,of,the,vector,potential,A||,=,b,·,A,where,b,=,B/,|B|,is,the,unit,vector,along,the,magnetic,field,B.,The,electrostatic,potential,is,calculated,from,a,fluid,vorticity,ω,,or,analogously,in,gyro-fluid,and,gyro-kinetic,models,,from,ion,and,gyro-centre,density,differ-,ences.,These,give,rise,to,an,equation,of,the,form,∇,·,(cid:16),n,B2,∇⊥φ,(cid:17),=,ω,(1),where,n,is,the,plasma,density,,which,varies,in,time,and,space.,This,system,is,often,simplified,by,replacing,n,by,a,constant.,By,analogy,with,buoyancy-driven,flows,,this,is,called,the,Boussinesq,approximation.,The,operator,∇⊥,=,∇−bb·∇,is,the,component,of,the,gradient,perpendicular,to,the,magnetic,field.,It,arises,because,the,ion,polarisation,drift,,which,ultimately,gives,rise,to,ω,,depends,on,2,the,electric,field,perpendicular,to,the,magnetic,field.,The,electric,field,parallel,to,the,magnetic,field,is,determined,by,quite,different,physics,,being,mainly,determined,by,the,electron,dynamics,rather,than,the,ions.,The,electromagnetic,potential,A||,=,b,·,A,is,related,to,the,current,along,the,magnetic,field:,B,=,∇,×,A,J,=,=,=,1,µ0,1,µ0,1,µ0,∇,×,B,∇,×,∇,×,A,(cid:2)∇,(∇,·,A),−,∇2A(cid:3),The,Coulomb,gauge,is,chosen,,setting,∇·A,=,0.,The,current,along,the,magnetic,field,is,therefore,J||,=,b,·,J,=,−,1,µ0,b,·,∇2A,(2),This,elliptic,operator,is,slightly,different,from,equation,1,above,,but,often,the,following,approximation,is,used:,b,·,∇2A,(cid:39),∇,·,(cid:0)∇A||,(cid:1),(3),In,addition,derivatives,of,A||,along,the,magnetic,field,are,often,neglected,,so,that,the,same,operator,as,equation,1,can,be,used.,The,potential,A||,appears,in,models,through,the,perturbed,magnetic,field,δB,=,∇,×,(cid:0)bA||,(cid:1),(4),It,also,appears,in,the,component,of,the,electric,field,parallel,to,the,magnetic,field:,E||,=,b,·,E,=,−b,·,∇φ,−,∂A||,∂t,(5),which,appears,in,an,Ohm’s,law,equation,for,the,current,along,the,magnetic,field.,3,Figure,1:,Domain,perpendicular,to,the,magnetic,field,in,a,torus.,For,illustration,only.,Left:,A,part,of,the,domain;,Right:,Space,filling,after,many,toroidal,turns.,2.1,Embedded,domain,The,elliptic,equation,to,be,solved,for,φ,in,equation,1,may,appear,to,be,a,2D,problem,,because,it,depends,on,the,component,of,the,electric,field,perpendicular,to,the,magnetic,field.,Indeed,if,the,magnetic,field,were,constant,in,space,then,this,would,be,a,2D,system.,In,a,tokamak,however,the,magnetic,field,varies,in,space,,and,generally,has,components,in,both,toroidal,and,poloidal,directions.,This,means,that,the,2D,domain,does,not,in,general,close,on,itself,,but,forms,a,spiral,which,fills,the,toroidal,volume.,This,is,illustrated,in,figure,1.,In,many,situations,the,tilt,of,the,2D,plane,in,the,toroidal,direction,can,be,neglected,,using,the,fact,that,variation,along,the,magnetic,field,(which,is,predominantly,toroidal),is,generally,small.,Some,important,exceptions,are:,•,Low-n,(toroidal,mode,number),instabilities,and,waves.,For,these,the,assumption,that,parallel,gradients,are,small,relative,to,perpendicular,can,break,down,,and,the,corrections,become,important.,•,In,low,aspect-ratio,(“spherical“),tokamaks,,the,pitch,of,the,magnetic,field,can,become,quite,large:,At,the,outboard,midplane,of,MAST,and,MAST-,U,,the,pitch,can,be,nearly,45o.,In,either,of,these,cases,solving,for,the,potential,may,require,a,3D,solve,,rather,than,a,decoupled,set,of,2D,solves.,It,seems,likely,that,in,most,cases,the,correc-,tions,should,be,small,,and,that,2D,solves,would,be,a,good,preconditioner,for,solving,the,full,3D,problem.,4,3,Tests,of,elliptic,solvers,As,for,many,components,of,scientific,high,performance,codes,,tests,need,to,address,both,the,correctness,and,performance,of,the,implementation.,3.1,Correctness,tests,In,all,correctness,tests,it,is,important,to,define,a,figure,of,merit.,This,should,be,quantitative,,and,sufficiently,well,defined,that,it,can,be,automated.,This,enables,a,pass/fail,criterion,to,be,established,,and,the,test,integrated,into,continuous,integration,workflows.,It,is,usually,useful,to,combine,both,a,measure,of,a,global,average,error,(such,as,root-mean-square,,l2,norm),,and,a,measure,which,is,sensitive,to,large,localised,errors,(such,as,maximum,absolute,error,,l∞,norm).,This,enables,both,the,overall,accuracy,of,the,scheme,to,be,assessed,,and,also,helps,identify,problems,with,specific,areas,of,the,mesh,such,as,boundaries.,Round-trip,tests:,The,simplest,tests,to,implement,are,those,which,use,a,forward,operator,to,check,the,accuracy,of,an,inversion,method.,This,type,of,test,is,often,fast,,and,useful,as,check,while,developing,the,code.,There,are,some,subtleties,in,this,which,can,lead,to,spurious,issues:,The,forward,operator,must,use,the,exact,same,discretisation,as,the,inverse,operator,being,tested.,The,forward,operator,itself,should,also,be,checked,for,correctness.,Manual,inspection,of,a,forward,operator,is,usually,more,straightforward,than,an,inverse,operator,,but,is,not,a,substitute,for,an,actual,quantitative,test.,Analytic,solution,tests:,In,simple,geometries,(such,as,slabs),,and,for,partic-,ular,choices,of,the,density,n,in,equation,1,,an,analytic,solution,can,be,found.,In,these,cases,the,numerical,solution,can,be,compared,against,the,analytic,,to,calculate,the,error.,To,properly,check,the,method,,a,convergence,test,should,be,performed,,using,multiple,grids,with,different,resolutions,,and,the,order,of,convergence,compared,against,the,theoretical,order,of,the,method,used.,Testing,the,order,of,convergence,in,this,way,can,be,challenging,in,some,cases:,Depend-,ing,on,the,system,,it,can,be,that,numerical,round-off,begins,to,affect,the,error,,before,the,error,enters,the,asymptotic,regime.,5,Manufactured,solution,tests:,Simple,analytic,tests,have,the,benefit,of,being,easy,to,implement,,but,are,typically,limited,in,their,thoroughness:,In,a,slab,,for,example,,many,metric,tensor,components,are,zero,which,in,a,realistic,simulation,would,be,non-zero.,In,that,case,the,slab,test,is,not,able,to,tell,whether,the,parts,of,the,code,which,depend,on,those,metric,components,are,implemented,correctly.,The,challenge,for,testing,is,that,realistic,cases,which,exercise,all,parts,of,the,code,typically,do,not,have,analytic,solutions,(otherwise,there,would,be,little,point,in,using,the,code).,Fortunately,since,the,system,to,be,tested,is,relatively,straightforward,,an,analytic,solution,can,be,chosen,(manufactured),for,φ,,differentiated,analytically,to,calculate,the,ω,function.,Values,for,ω,can,then,be,calculated,on,a,grid,,inverted,and,compared,to,the,original,expression.,As,part,of,this,,a,non-trivial,geometry,needs,to,be,generated,analytically.,This,need,not,necessarily,be,of,the,same,form,as,real,simulations,(e.g.,a,tokamak),,provided,that,it,exercises,the,same,parts,of,the,code,(preferably,,all,the,code).,3.2,Performance,tests,The,elliptic,solver,is,often,the,main,bottleneck,to,parallel,scaling,of,plasma,physics,codes,,certainly,in,existing,codes,such,as,BOUT++.,This,is,because,it,involves,global,communications,,or,at,least,communications,over,a,large,region,of,the,domain,,and,this,inversion,is,done,frequently,as,part,of,the,time,advance.,Elliptic,solvers,are,critical,to,many,areas,of,scientific,computing,,certainly,not,unique,to,plasma,physics.,Experience,with,BOUT++,however,,has,been,that,the,plasma,use,case,is,different,from,others,in,ways,which,are,important,for,performance:,Libraries,are,often,optimised,for,solving,very,large,systems,of,equations,(millions,of,d.o.f),,for,applications,where,a,relatively,small,number,of,such,large,systems,may,need,to,be,solved.,In,typical,plasma,turbulence,simulations,,however,,the,size,of,each,system,may,be,relatively,small,(104,d.o.f,for,a,2D,system;,100s,of,degrees,of,freedom,per,solve,if,decomposed,into,1D,systems),,but,a,turbulence,simulation,may,require,105,or,106,such,solves.,Since,the,systems,to,be,solved,are,part,of,the,time,evolution,,these,105,or,106,solves,are,essentially,serial,unless,parallel-in-time,techniques,are,used,,and,must,be,solved,efficiently.,Since,the,quantity,being,solved,for,is,evolving,in,time,,solutions,tend,to,be,close,to,previous,solutions,,and,so,iterative,schemes,tend,to,be,effective.,6,It,is,more,important,that,performance,tests,use,problems,which,are,close,to,those,which,will,be,solved,in,production,systems,,than,it,is,for,correctness,tests.,Different,problem,sizes,hit,different,thresholds,in,cache,sizes,,message,sizes,,work,intensity,,network,capacity,etc.,Problems,in,slab,geometry,,which,are,simpler,to,solve,than,realistic,situations,,may,have,different,convergence,characteristics,,and,not,exercise,the,preconditioner,,changing,the,proportion,of,time,spent,in,parts,of,the,solver.,As,discussed,above,,a,characteristic,of,typical,problems,is,that,they,are,solving,systems,which,start,from,a,generally,good,initial,guess,,the,value,at,the,previous,time,step.,This,should,be,taken,into,account,in,the,performance,testing.,Specific,solution,algorithms,,their,scaling,,theoretical,and,measured,perfor-,mance,are,not,discussed,in,this,report,beyond,the,comments,above.,These,will,be,addressed,in,task,1.2,and,1.3.,4,Tests,in,slab,geometries,A,series,of,slab,simulations,of,gradually,increasing,complexity,can,be,used,in,build,a,code,up,step-by-step,,and,help,to,identify,problems,early,before,moving,on,to,more,complex,systems.,•,A,2D,domain,,with,fixed,(Dirichlet),boundary,conditions,on,all,sides,,and,the,magnetic,field,directed,perpendicular,to,the,domain.,This,is,a,simple,problem,which,is,probably,included,as,an,example,in,most,finite,difference,or,finite,element,packages.,•,The,magnetic,field,vector,b,is,not,perpendicular,to,the,plane,the,potential,φ,is,being,solved,on.,This,introduces,some,non-zero,metric,terms,,and,can,be,made,more,challenging,by,varying,the,angle,of,the,magnetic,field,with,position.,•,The,boundaries,can,be,modified,so,that,the,mesh,is,periodic,in,one,di-,rection,,but,continues,to,have,boundaries,in,the,other.,In,a,tokamak,the,periodic,direction,would,correspond,to,the,poloidal,or,toroidal,directions,(depending,on,the,coordinate,system,chosen),,and,the,other,direction,to,the,radial,direction,from,inside,(core),to,outside,(edge).,7,•,A,variation,on,the,above,is,to,introduce,a,corner,into,the,problem,,so,that,the,domain,is,periodic,over,part,of,the,boundary,(in,the,core),,but,has,fixed,boundaries,over,the,rest,of,the,boundary,(the,scrape-off,layer).,4.1,Slab,with,a,pole,It,was,found,in,developing,a,BOUT++,implementation,,that,the,above,slab,ge-,ometries,were,not,sufficiently,challenging,to,be,able,to,exercise,a,GPU,precon-,ditioner,(Hypre).,To,make,the,problem,more,challenging,,and,closer,to,realistic,tokamak,geometry,,a,test,case,was,developed,which,has,a,pole,in,the,coordinate,system,close,to,(but,outside),the,mesh,edge.,This,mimics,the,singularity,at,the,X-point,which,occurs,in,field-aligned,coordinates.,In,field-aligned,coordinates,one,of,the,coordinates,is,aligned,with,the,magnetic,field.,In,this,case,the,y,coordinate,basis,vector,ey,is,aligned,to,B,,so,therefore,the,gradients,of,the,x,and,z,coordinates,are,perpendicular,to,B.,A,choice,used,in,BOUT++,and,other,plasma,simulation,codes,is,a,Clebsch,coordinate,system,,characterised,by,B,=,∇z,×,∇x.,For,details,see,the,BOUT++,manual,section,on,field-aligned,coordinates,[1].,The,radius,from,a,pole,,r,is,used,to,set,non-zero,metric,tensor,components,gxx,=,ex,·,ex,=,1/r2,gyy,=,ey,·,ey,=,1,+,1/r2,gzz,=,ez,·,ez,=,1,gyz,=,ey,·,ez,=,1/r,(6),(7),(8),(9),and,so,coordinate,Jacobian,J,=,1/r.,In,BOUT++,the,input,file,for,this,configuration,is,shown,in,figure,2.,5,Tests,in,tokamak,geometries,These,tests,approach,realistic,production,cases,,introducing,a,toroidal,geometry,with,sheared,magnetic,field,in,doubly-periodic,domains,,and,then,by,introducing,8,1,#,Mesh,size,2,Lx,=,10,3,Ly,=,10,4,5,#,Number,of,grid,cells,6,nx,=,20,7,ny,=,32,8,9,#,mesh,spacing,10,dx,=,Lx,/,(,nx,-,4),#,Account,for,4,guard,cells,in,X,11,dy,=,Ly,/,ny,12,13,#,Location,of,the,pole,in,the,coordinates,14,pole_x,=,Lx,+,1.0,15,pole_y,=,Ly,+,1.0,16,17,#,Distance,from,the,pole,18,#,Note,here,",x,",is,normalised,to,[0,,1],on,the,grid,;,",y,",normalised,to,[0,,2*,pi,],19,r,=,sqrt,((,pole_x,-,x,*,Lx,),^2,+,(,pole_y,-,y,*,Ly,/,(2*,pi,),),^2),20,21,#,This,mimics,the,metric,tensor,close,to,the,X,-,point,in,a,tokamak,22,#,by,here,setting,poloidal,field,Bp,~,r,23,24,g11,=,r,^2,25,g22,=,1,26,g33,=,1,+,1,/,r,^2,27,g12,=,0,28,g13,=,0,29,g23,=,-1/,r,30,31,J,=,1,/,r,32,33,g_11,=,1,/,r,^2,34,g_22,=,1,+,1,/,r,^2,35,g_33,=,1,36,g_12,=,0,37,g_13,=,0,38,g_23,=,1,/,r,Figure,2:,Part,of,a,BOUT++,input,file,for,a,slab,with,a,pole.,See,BOUT++,manual,for,details,of,the,code,[2].,an,X-point,into,the,domain.,For,each,geometry,two,kinds,of,tests,can,be,performed:,•,A,single,solve,(for,correctness),,or,small,number,of,solves,with,slowly,varying,input,(for,performance).,•,Time-evolving,a,relatively,simple,system,of,equations,,in,which,the,elliptic,solve,is,a,key,part.,The,purpose,of,this,is,to,identify,potential,issues,with,9,slowly,accumulating,errors.,These,errors,may,be,below,tolerance,in,a,single,solve,,but,interact,with,the,time,evolving,system,in,a,way,which,leads,to,numerical,instability.,5.1,Tokamak,geometries,Tokamak,grids,are,typically,generated,from,a,numerical,equilibrium,,and,typ-,ically,those,are,provided,from,experiment,at,low,resolution,(e.g.,64x64,for,the,whole,poloidal,cross-section).,Higher,resolution,inputs,can,be,generated,using,a,free,boundary,Grad-Shafranov,solver,,but,generating,a,sequence,of,sim-,ulation,meshes,from,Grad-Shafranov,solutions,for,a,convergence,test,remains,challenging.,Not,conceptually,difficult,,but,the,capability,to,do,it,has,not,been,implemented,and,would,require,some,considerable,effort,to,do,in,a,rigorous,way.,For,convergence,tests,,the,easiest,approach,would,seem,to,be,to,use,an,analytic,solution,,either,to,the,Grad-Shafranov,solution,itself,,or,a,simplified,solution,which,is,not,a,Grad-Shafranov,solution,but,shares,important,characteristics,(such,as,an,X-point).,5.1.1,Numerical,tokamak,equilibria,If,an,analytic,equilibrium,and,convergence,to,small,tolerances,is,not,required,,then,numerical,solutions,are,available,for,many,different,tokamaks,,both,real,and,conceptual.,A,common,starting,point,is,circular,cross-section,geometries,,which,don’t,have,an,X-point.,Examples,include,the,‘cbm18’,series,of,equilibria,generated,by,P.Snyder,(GA).,Alternatively,an,equilibrium,can,be,generated,by,a,free-boundary,equilibrium,code,such,as,EFIT,,EFIT++,,Helena,,or,FreeGS[3].,5.1.2,Analytic,tokamak,equilibria,The,Grad-Shafranov,equation,is,a,nonlinear,partial,differential,equation,,and,so,finding,analytic,solutions,is,non-trivial.,Some,have,however,been,found:,The,Soloviev,solutions,are,widely,used,,but,only,include,closed,flux,surfaces,(the,plasma,core),,not,the,separatrix,and,scrape-off,layer.,In,addition,these,10,solutions,typically,have,jumps,or,current,sheets,at,the,edge,which,make,them,problematic,to,extend,into,the,vacuum,Cerfon,&,Freidberg,found,a,set,of,analytic,equilibria,which,include,a,separatrix,and,vacuum,region,outside,the,plasma,[4].,That,algorithm,has,been,imple-,mented,in,Python,by,John,Omotani,[5].,This,could,be,used,as,input,to,a,mesh,generator,for,convergence,studies.,6,Time-evolving,systems,These,are,simple,systems,of,equations,which,are,intended,to,be,relatively,quick,to,simulate,,but,which,can,help,identify,problems,with,numerical,methods,or,boundary,conditions,at,an,early,stage.,Since,we,wish,to,identify,numerical,instabilities,and,error,accumulation,,we,choose,systems,which,are,stable:,These,systems,of,equations,support,waves,which,are,either,stable,or,damped.,Any,growing,mode,can,therefore,be,identified,as,a,numerical,artefact.,6.1,Shear,Alfv`en,wave,This,is,an,electromagnetic,wave,along,the,magnetic,field.,Due,to,the,variation,of,the,magnetic,field,in,a,tokamak,,there,is,a,radial,(across,the,field),phase,velocity,,and,an,initially,coherent,mode,will,tend,to,mix,and,dissipate,due,to,the,model,and,numerical,damping.,The,set,of,equations,evolves,the,vorticity,U,and,parallel,vector,potential,A||.,It,appears,in,normalised,form,as:,∂U,∂t,∂A||,∂t,=,∇||j||,=,−∂||φ,−,ηj||,+,µ||e∇2,||j||,2,βe,∇2,⊥A||,j||,=,−,∇2,⊥φ,=,U,11,(10),(11),(12),(13),where,∇||f,=,∇,·,(bf,),is,the,divergence,of,a,flow,along,the,magnetic,field,,and,∂||,=,b,·,∇,is,the,gradient,along,the,magnetic,field.,The,factor,βe,is,the,ratio,of,electron,pressure,to,magnetic,pressure,,and,is,typically,around,10−4,in,the,plasma,edge.,Dissipation,is,included,in,the,form,of,resistivity,η,and,parallel,electron,viscosity,µ||e.,The,resistivity,would,typically,be,small,(<,10−2),,and,electron,viscosity,much,smaller,than,that,(usually,neglected).,6.1.1,Energy,conservation,The,equations,for,this,wave,contain,only,energy,sinks,(dissipation),,and,so,oscillations,can,only,grow,if,there,are,sources,of,energy,from,either,numerical,instability,or,boundary,fluxes.,To,calculate,the,energy,in,the,system,,multiply,the,vorticity,equation,by,φ,and,the,A||,equation,by,j||,,then,integrate,over,the,simulation,volume.,φU,=,φ∇,·,∇⊥φ,=,∇,·,(φ∇⊥φ),−,∇φ,·,∇⊥φ,(cid:125),(cid:124),(cid:123)(cid:122),|∇⊥φ|2,∂,∂t,(φU,),=,φ,+,U,∂U,∂φ,∂t,∂t,=,φ∇||j||,+,∇2,∇2,⊥φ,∂φ,∂t,=,∇,·,(cid:18),∂φ,∂t,∂φ,∂t,⊥φ,(cid:19),∇⊥φ,−,∇,∂φ,∂t,·,∇⊥φ,(cid:125),(cid:123)(cid:122),∂t,(|∇⊥φ|2),∂,(cid:124),1,2,(14),(15),(16),(17),and,so,(cid:104),∇,·,(φ∇⊥φ),−,|∇⊥φ|2(cid:105),=,φ∇||j||,+,∇,·,∂,∂t,(cid:18),∂φ,∂t,(cid:19),∇⊥φ,−,1,2,∂,∂t,(cid:16),|∇⊥φ|2(cid:17),(18),Putting,these,together,gives,an,equation,for,the,evolution,of,the,kinetic,energy,in,the,E×B,motion:,1,2,∂,∂t,(cid:16),|∇⊥φ|2(cid:17),=,−φ∇||j||,+,∇,·,φ∇⊥,(cid:18),(cid:19),∂φ,∂t,(19),On,the,left,is,the,energy,in,the,E×B,motion,,since,the,factors,of,ion,mass,and,density,are,constant,here,,and,have,been,normalised,out.,The,first,term,on,the,right,is,a,transfer,of,energy,from,electromagnetic,energy,,and,the,second,term,is,12,a,flux,of,energy,from,the,boundary.,Integrating,this,equation,over,a,volume,,the,divergence,term,on,the,right,will,become,a,surface,integral,over,the,boundary:,∂,∂t,(cid:90),V,1,2,|∇⊥φ|2,dV,=,−,(cid:90),V,φ∇||j||dV,+,(cid:73),S,φ∇⊥,∂φ,∂t,·,dS,(20),and,so,the,flux,of,energy,through,the,boundary,is,zero,if,φ,is,zero,,or,if,the,component,of,∇⊥φ,normal,to,the,boundary,is,constant,in,time.,Similarly,,the,A||,equation,gives:,j||A||,=,2,βe,A||∇2,⊥A||,=,2,βe,∇,·,(cid:0)A||∇⊥A||,(cid:1),−,2,βe,∇A||,·,∇⊥A||,(cid:125),(cid:123)(cid:122),(cid:124),|∇⊥A|||2,∂,∂t,(cid:0)j||A||,(cid:1),=,j||,=,j||,+,A||,∂A||,∂t,(cid:0)∂||φ,+,ηj||,∂j||,∂t,(cid:1),+,(cid:18),A||∇,·,∇⊥,2,βe,(cid:19),∂A||,∂t,(cid:19),=,j||∂||φ,+,ηj2,||,+,(cid:18),∇,·,A||∇⊥,2,βe,∂A||,∂t,−,2,βe,1,2,∂,∂t,(cid:12),(cid:12)∇⊥A||,(cid:12),(cid:12),2,(24),(21),(22),(23),and,so:,∂,∂t,(cid:20),2,βe,∇,·,(cid:0)A||∇⊥A||,(cid:1),−,2(cid:21),(cid:12),(cid:12)∇⊥A||,(cid:12),(cid:12),2,βe,=,j||∂||φ+ηj2,||+,(cid:18),∇·,A||∇⊥,∂A||,∂t,2,βe,(cid:19),−,2,1,2,βe,(25),∂,∂t,(cid:12),(cid:12)∇⊥A||,(cid:12),(cid:12),2,1,βe,∂,∂t,(cid:12),(cid:12)∇⊥A||,(cid:12),(cid:12),2,=,−j||∂||φ,−,ηj2,||,+,2,βe,∇,·,(cid:18),∂A||,∂t,(cid:19),∇⊥A||,(26),Integrating,over,the,volume:,∂,∂t,(cid:90),V,1,βe,(cid:12),(cid:12)∇⊥A||,(cid:12),(cid:12),2,dV,=,−,(cid:90),V,j||∂||φdV,−,(cid:90),V,ηj2,||dV,+,(cid:73),S,2,βe,∂A||,∂t,∇⊥A||,·,dS(27),Like,the,E×B,energy,equation,(eq,20),,the,flux,of,energy,through,the,boundary,is,zero,if,A||,=,const,or,∇⊥A||,·,S,=,0.,The,exchange,of,energy,between,E×B,and,electromagnetic,forms,is,through,13,φ∇||j||,and,j||∂||φ.,These,should,balance,for,energy,to,be,conserved:,φ∇||j||,=,φ∇,·,(cid:0)bj||,(cid:1),=,∇,·,(cid:0)bφj||,(cid:1),−,j||,b,·,∇φ,(cid:124),(cid:123)(cid:122),(cid:125),∂||φ,(28),Hence,there,are,no,fluxes,of,energy,through,the,parallel,boundary,if,j||,or,φ,are,zero,at,the,boundary.,From,this,we,conclude:,•,φ,=,0,or,∇⊥φ,·,S,=,const,for,conservation,of,energy,•,A||,=,const,or,∇⊥A||,·,S,=,0,for,conservation,of,energy,•,There,is,no,requirement,that,the,form,of,∇2,⊥,in,vorticity,and,j||,equations,is,the,same,,since,the,only,term,which,is,common,to,both,equations,is,φ∇||j||,↔,j||∂||φ,6.2,Geodesic,Acoustic,Wave,The,Geodesic,Acoustic,Mode,(GAM),is,an,axisymmetric,(n,=,0),sound,wave,which,occurs,through,oscillations,in,the,radial,current.,It,involves,an,up-down,(m,=,1),asymmetry,in,the,density,,together,with,a,potential,which,is,approx-,imately,constant,on,flux,surfaces,(m,=,0).,It,can,arise,from,a,simple,system,of,equations,,but,because,it,involves,currents,across,the,magnetic,field,,it,exer-,cises,the,radial,boundary,conditions.,If,an,annulus,is,simulated,,so,that,there,is,a,boundary,to,the,inner,core,,then,treatment,of,that,core,boundary,can,be,important,,to,ensure,that,there,is,no,net,current,,or,that,radial,boundary,layers,don’t,form.,6.2.1,Simplified,model,The,GAM,oscillation,arises,because,radial,gradients,of,zonal,(n,=,0,,m,(cid:39),0),electrostatic,potential,causes,E×B,advection,of,density,in,the,poloidal,direction.,Because,the,magnetic,field,strength,varies,between,inboard,and,outboard,sides,,E×B,advection,is,faster,on,the,outboard,side,than,the,inboard.,This,leads,to,rarefaction,and,compression,at,the,top,and,bottom,of,the,device,(depending,14,on,flow,direction).,The,change,in,pressure,at,top,and,bottom,of,the,device,alters,the,diamagnetic,current,across,the,magnetic,field,,which,then,modifies,the,electrostatic,potential.,We,require,equations,for,current,continuity,and,density,n:,∇,·,(cid:20),min,B2,∂∇⊥φ,∂t,(cid:21),(cid:20),=,∇,·,p∇,×,(cid:19)(cid:21),(cid:18),b,B,+,∇||j||,∂n,∂t,(cid:20),=,−∇,·,n,(cid:21),b,×,∇φ,B,(29),(30),The,first,(vorticity),equation,includes,a,current,along,the,magnetic,field,,j||.,This,is,not,required,to,derive,the,wave,analytically,,but,is,required,in,a,numerical,simulation,to,ensure,that,φ,remains,approximately,constant,along,the,magnetic,field.,For,this,purpose,a,simple,electrostatic,Ohm’s,law,can,be,used:,j||,=,−,1,η,b,·,∇φ,(31),Decreasing,η,will,make,the,potential,relax,more,quickly,to,a,constant,along,flux,surfaces.,The,usual,approach,to,solving,this,is,to,use,∇,×,b,B,(cid:39),2,B,b,×,κ,and,to,split,the,E,×,B,advection,term,into,a,divergence-free,advection,term,,and,a,divergence,term:,∂n,∂t,=,−,1,B,b,×,∇φ,·,∇n,−,n∇,·,(cid:19),b,×,∇φ,(cid:18),1,B,then,approximate,∇,·,(cid:18),1,B,(cid:19),b,×,∇φ,(cid:39),2,B,b,×,κ,·,∇φ,(32),(33),It,is,usual,to,neglect,the,poloidal,derivative,(y),terms,in,the,E,×,B,advection,operator.,In,Clebsch,coordinates,this,term,looks,like,1,B,b,×,∇φ,·,∇n,=,∂φ,∂x,∂φ,∂z,−,∂φ,∂z,∂φ,∂x,15,For,axisymmetric,flows,the,z,derivatives,are,zero,,so,this,term,vanishes,,leaving,only,the,compression,term.,This,leaves,a,minimal,GAM,model:,∂Ω,∂t,∂n,∂t,=,2,B,b,×,κ,·,∇p,+,∇||j||,=,−n,2,B,min0,B2,0,Ω,=,b,×,κ,·,∇φ,∇2,⊥φ,(34),(35),(36),where,n0,and,B0,are,constants,in,space,and,time,,and,an,isothermal,approxi-,mation,is,used,here:,p,=,enT0,(37),Note,also,that,here,φ,is,assumed,to,be,approximately,constant,on,flux,surfaces,,which,will,need,to,be,enforced,numerically,using,an,Ohm’s,law.,6.2.2,Analytic,solution,Starting,with,the,density,equation,,we,look,for,solutions,of,the,form,and,potential,φ:,n,(x,,θ,,t),=,ˆn,(θ),eikx−iwt,φ,(x,,θ,,t),=,ˆφeikx−iwt,Linearising,equation,35,,assuming,a,simple,circular,cross-section,,large,aspect-,ratio,,and,keeping,only,the,poloidal,flow,we,get,b,×,κ,·,∇,→,1,R,sin,θ,∂,∂x,ˆn,=,n0,2k,BRω,sin,θ,ˆφ,(38),(39),Hence,if,we,assume,ˆφ,is,independent,of,θ,then,ˆn,has,a,sin,θ,dependence.,The,parallel,current,term,will,act,to,equalise,potential,over,a,flux,surface,,but,provided,this,occurs,sufficiently,rapidly,we,can,remove,it,from,the,analysis,and,assume,that,φ,is,constant,on,flux,surfaces.,To,remove,the,parallel,current,term,16,from,the,vorticity,equation,,average,over,a,flux,surface,by,defining,(cid:104)·(cid:105),so,that,(cid:10)∇||,(cid:11),=,0.,For,large,aspect,ratio:,(cid:73),·dθ,(cid:104)·(cid:105),=,This,gives:,−iω,min0,B2,(cid:104)Ω(cid:105),=,∂,∂t,(cid:69),(cid:0)−k2(cid:1),(cid:68),ˆφ,ωk2,min0,B2,(cid:28),(cid:29),b,×,κ,·,∇p,(cid:28),2,B,eT0,BR,k2en0T0,ωB2R2,sin2,θ,sin,θ,(ik),ˆn,(cid:28),2,4,(cid:29),(cid:29),(40),(41),(42),=,=,Most,terms,cancel,,leaving,eT0,mi,i.e.,the,frequency,depends,on,the,sound,speed,cs,=,(cid:112)eT0/mi,and,major,radius,R.,This,is,the,correct,dependency,for,GAMs,in,the,simple,electrostatic,,large,2,R2,ω2,=,aspect,ratio,limit,when,parallel,flows,are,neglected.,6.2.3,Mass,and,energy,conservation,For,accurate,and,stable,numerical,simulations,the,mass,and,energy,of,the,system,should,be,conserved,over,long,times.,Starting,by,multiplying,the,vorticity,equation,by,φ,φΩ,=,φ,min0,B2,0,∇2,⊥φ,=,min0,B2,0,(cid:104),∇,·,(φ∇⊥φ),−,|∇⊥φ|2(cid:105),∂,∂t,(φΩ),=,φ,=,φ,∂Ω,∂t,2,B,+,Ω,∂φ,∂t,b,×,κ,·,∇p,+,φ∇||j||,+,Ω,(43),(44),∂φ,∂t,Ω,∂φ,∂t,=,min0,B2,0,(cid:20),∇,·,(cid:18),∂φ,∂t,(cid:19),∇⊥φ,−,1,2,∂,∂t,(cid:21),|∇⊥φ|2,17,∂,∂t,(cid:20),1,2,min0,B2,0,(cid:21),|∇⊥φ|2,=,−φ,2,B,min0,B2,0,+,b,×,κ,·,∇p,−,φ∇||j||,(45),(cid:20),∂,∂t,∇,·,(φ∇⊥φ),−,∇,·,(cid:18),∂φ,∂t,(cid:19)(cid:21),∇⊥φ,(46),This,equation,describes,the,change,of,E,×,B,kinetic,energy.,The,first,line,(eq,45),contains,transfer,terms,from,other,forms,of,energy,,whilst,the,second,term,describes,fluxes,from,the,boundaries.,By,setting,either,φ,=,0,or,∇⊥φ·S,=,0,at,the,boundary,the,fluxes,go,to,zero,at,the,boundary.,Many,choices,for,the,parallel,current,are,possible,,but,a,simple,form,is,E||,=,−∂||φ,=,ηj||,where,η,is,the,resistivity.,In,this,case,∇||,(cid:0)φj||,(cid:1),=,φ∇||j||,+,j||∂||φ,=,φ∇||j||,−,ηj2,||,−φ∇||j||,=,−ηj2,||,−,∇||,(cid:0)φj||,(cid:1),(47),(48),(49),(50),and,so,this,term,is,always,a,sink,of,energy,,and,boundary,fluxes,go,to,zero,if,φ,=,0,or,j||,=,0,at,the,boundary.,The,curvature,transfer,term,can,be,derived,by,starting,from,(cid:18),∇,·,φp,2,B,(cid:19),b,×,κ,=,φ,2,B,b,×,κ,·,∇p,+,p,2,B,b,×,κ,·,∇φ,+,pφ∇,·,(cid:19),b,×,κ,(51),(cid:18),2,B,∂,∂t,(cid:20),1,2,min0,B2,0,(cid:21),|∇⊥φ|2,=,p,b,×,κ,·,∇φ,−,ηj2,||,(52),∇,·,(φ∇⊥φ),−,∇,·,(cid:18),∂φ,∂t,(cid:19)(cid:21),∇⊥φ,(53),−,∇||,(cid:0)φj||,(cid:1),(cid:19),b,×,κ,(cid:19),b,×,κ,(54),(55),2,B,min0,B2,0,(cid:18),+,−,∇,·,φp,+,pφ∇,·,(cid:20),∂,∂t,2,B,(cid:18),2,B,18,Multiplying,equation,35,by,the,constant,eT0,we,get,∂p,∂t,=,−p,2,B,b,×,κ,·,∇φ,(56),It,can,be,seen,that,this,term,balances,the,first,term,on,the,right,of,equation,52.,The,total,energy:,E,=,(cid:90),V,dV,(cid:20),1,2,min0,B2,0,(cid:21),|∇⊥φ|2,+,p,(57),is,conserved,(apart,from,resistive,losses),so,long,as,the,terms,in,equations,53,,54,and,55,vanish.,1.,Equations,53,and,54,are,divergences,,and,all,go,to,zero,at,the,boundaries,if,φ,=,0,at,the,boundary.,2.,Equation,55,is,not,a,divergence,,so,can,produce,sources,and,sinks,of,energy,in,the,domain,,not,just,at,the,boundary.,Since,p,and,φ,can,be,arbitrary,functions,,the,curvature,vector,must,satisfy,∇,·,(cid:18),2,B,(cid:19),b,×,κ,=,0,(58),The,toroidal,component,of,this,curvature,vector,2,B,b,×,κ,doesn’t,affect,the,conservation,properties;,the,radial,x,component,is,essential,for,the,GAM.,Either,the,x,component,has,to,be,constant,,or,the,y,component,must,also,be,included.,7,Conclusions,A,set,of,test,cases,have,been,described,,starting,from,simple,slabs,and,pro-,gressing,first,to,more,complex,geometries,,and,then,to,integrating,the,elliptic,solver,as,part,of,a,time-evolving,system.,Criteria,for,correctness,have,been,specified:,l2,and,l∞,error,norms,for,tests,with,an,analytic,solution,(manufac-,tured,in,complex,cases).,For,time-evolving,problems,energy,conservation,is,a,useful,measure,,because,it,has,a,strong,connection,to,numerical,stability,and,physical,correctness,of,the,result.,These,tests,will,be,used,in,future,reports,to,test,elliptic,solver,implementations.,19,8,References,[1],BOUT++,contributors.,BOUT++,manual:,Field-aligned,coor-,dinates.,coordinates.html#id1.,https://bout-dev.readthedocs.io/en/latest/user_docs/,[2],BOUT++,contributors.,BOUT++,manual.,https://bout-dev.,readthedocs.io/.,[3],FreeGS,contributors.,FreeGS:,Free,boundary,Grad-Shafranov,solver,in,python.,https://github.com/bendudson/freegs.,[4],Antoine,J.,Cerfon,and,Jeffrey,P.,Freidberg.,“One,size,fits,all”,analytic,solutions,to,the,Grad–Shafranov,equation.,Phys.,Plasmas,,17:032502,,2010,,doi:10.1063/1.3328818.,[5],John,Omotani.,Cerfon-freidberg,geometry,generator,in,python.,https:,//github.com/johnomotani/CerfonFreidbergGeometry.,20

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