APARTICLE-IN-CELLMETHODFORTHESIMULATIONOFPLASMASBASEDON ANUNCONDITIONALLYSTABLEWAVEEQUATIONSOLVER By EricMatthewWolf ADISSERTATION Submittedto MichiganStateUniversity inpartialentoftherequirements forthedegreeof AppliedMathematics-DoctorofPhilosophy 2015 ABSTRACT APARTICLE-IN-CELLMETHODFORTHESIMULATIONOFPLASMAS BASEDONANUNCONDITIONALLYSTABLEWAVEEQUATION SOLVER By EricMatthewWolf Inthisdissertation,wepresentaparticle-in-cellmethodforthesimulationofplasmas basedonanunconditionallystablesolverforthesecond-orderscalarwaveequation,that is,awaveequationsolverthatisnotsubjecttoaCourant-Friedrichs-Lewy(CFL)stability restriction,typicalofexplicitmethods,whilemaintainingacomputationalcostandcode complexitycomparabletosuchexplicitsolvers.Thispermitstheuseofatimestepsize manytimeslargerthanallowedbywidely-usedexplicitmethods. ACKNOWLEDGMENTS Thisworkwasdoneinanextremelysupportiveatmosphere,forwhichIamverygrateful. Ithankmyfamilyforalloftheirloveandsupportthroughoutmylifeandeducation.I thankmyadvisor,Prof.AndrewChristlieb,forhispatientguidanceovermanyyears, andforbuildinghisresearchgroupasanencouragingenvironment,inwhichthereareso manyknowledgableandfriendlypostdocsandgraduatestudentsdoingexcitingresearch inmanytareas.Ithankmyfellowgraduatestudentsandthepostdocsforallof yourinteractionandhelp.IespeciallythankProf.MatthewCausley,whoinhistimeasa postdocinourgroupderivedmuchofthemathematicsunderlyingtheimplicitwavesolver motivatingthisdissertation,andwhowasaninvaluablementortomeandothersinthe group.IalsothankDr.AndrewGreenwood,whomentoredmeduringasummerinternship andtaughtmeaboutparticle-in-cellmethods,andDr.MatthewBettencourt,whomentored meinanothersummerinternshipandgavemevaluableguidanceasIstarteddevelopingthe particle-in-cellmethodpresentedinthisdissertation. AndIthankmygirlfriendKazuko,mybestfriendformanyyears,whowhilethiswork wasdoneenduredwithmelongyearsofseparation.Shegavemeherpatience,love,and support. iii TABLEOFCONTENTS LISTOFTABLES .................................... vi LISTOFFIGURES ................................... vii Chapter1Introduction ............................... 1 1.1Overview......................................4 1.2MathematicalModels...............................6 1.3NondimensionalizationandAsymptoticAnalysis................8 1.3.1ElectrostaticCase.............................9 1.3.2ElectromagneticCase...........................12 Chapter2Descriptionofthewaveequationsolver .............. 14 2.1OriginandMitigationoftheCFLRestriction.................14 2.2DerivationofSecond-OrderSemi-DiscreteSchemesWiththeMethodofLines Transpose.....................................17 2.2.1DispersiveSemi-DiscreteScheme....................18 2.2.2eSemi-DiscreteScheme.....................18 2.3Semi-DiscreteVonNeumannAnalysisoftheSecondOrderSchemes.....19 2.3.1DispersiveScheme............................20 2.3.2eScheme.............................20 2.4SolutionoftheMoHelmholtzEquation..................22 2.4.1DimensionalSplitting...........................23 2.4.2IntegralSolutionMethodin1D.....................24 2.4.3FastNumericalEvaluationofthe1DConvolutionOperator......25 2.5FullyDiscreteVonNeumannAnalysis......................27 2.6DissipationinCenteredSchemes...................33 2.6.1Dissipationin1D........................33 2.6.2Dissipationin2D........................35 2.7DerivationofaFourth-OrderDissipativeNewmarkScheme..........36 2.8EmbbededBoundaryMethodforNeumannBoundaryConditions......40 2.8.1DescriptionoftheTwo-PointBoundaryCorrectionMethod......42 2.8.2DescriptionoftheEmbeddedBoundaryMethodandProofofConver- genceoftheIterativeSolutionin1D..................42 2.8.3DescriptionoftheMethodin2D....................48 2.9NumericalExamples...............................54 2.9.1RectangularCavityWithDomainDecomposition...........54 2.9.2DoubleCircleCavity...........................56 2.9.3PeriodicSlitGrating.....................58 2.9.4BesselModewithNeumannBoundaryConditions...........60 iv Chapter3Descriptionoftheparticle-in-cellmethod ............. 63 3.1ParticleWeightingScheme............................63 3.2MethodforControllingDivergenceError....................64 3.3ParticleEquationsofMotion...........................65 3.4FastConvolutionAlgorithmForParticleSources................66 3.4.1FastConvolutionAlgorithmin1D....................66 3.4.1.1LocalDepositStepin1D...................67 3.4.1.2RecursiveSweepStepin1D..................69 3.4.2FastConvolutionAlgorithmin2D....................70 3.4.2.1LocalDepositStepin2D...................70 3.4.2.2RecursiveSweepStepin2D..................72 3.5ParticleBoundaryConditions..........................73 3.5.11DPeriodicBC..............................74 3.5.22DPeriodicBC..............................74 3.5.31DDirichletBC..............................74 3.5.3.1AnalysisoftheFirstApproach................75 3.5.3.2AnalysisoftheSecondApproach...............77 3.5.3.3ComparisonoftheApproaches................78 3.5.42DDirichletBC..............................80 3.6NumericalResults.................................80 3.6.1ElectrostaticTestProblems.......................80 3.6.1.1ColdPlasmaLangmuirWave.................81 3.6.1.2TwoStreamInstability.....................81 3.6.1.3LandauDamping........................85 3.6.1.4SheathFormationinaBounded1DPlasma.........85 3.6.2ElectromagneticTestProblems.....................86 3.6.2.1BennettPinchProblem....................86 3.6.2.2MardahlBeamProblem....................91 Chapter4Conclusion ................................ 101 BIBLIOGRAPHY .................................... 102 v LISTOFTABLES Table2.1:Tablesummarizingthemainalgorithmicaspectsoftheimplicitwave solver...................................53 Table2.2:tstudyforrectangularcavitywithDirichletBCusingdo- maindecomposition.Here, c =1, m = n =0,and L x = L y =1...55 Table2.3:tstudyforrectangularcavitywithNeumannBCusing domaindecomposition.Here, c =1, m = n =0,and L x = L y =1..56 Table2.4:tstudyforthedoublecirclecavitywithDirichletBC. Forthenumericalreferencesolution, x =2 : 1875 10 4 , y = 1 : 3542 10 4 ,and t =2 : 7083 10 4 .................58 vi LISTOFFIGURES Figure2.1:Boundarygeometryin1D........................44 Figure2.2:Boundarygeometryin2D........................55 Figure2.3:Rectangularcavitywithdomaindecomposition............55 Figure2.4:Doublecirclegeometry..........................56 Figure2.5:Evolutionofthedoublecirclecavityproblem..............57 Figure2.6:Periodicslitgratinggeometry...............59 Figure2.7:Evolutionoftheslitgratingproblem,withaperturewidth a =0 : 1,gratingperiodicity L y = d =1,andwavespeed c =1.The CFLisat2.............................61 Figure2.8:Anexampleoftheembeddedboundarygridused.Theredcircledex- teriorgridpointsaretheendpointswhereavalueistobecalculated viatheinterpolationprocedure.Theredcrossesarethepointswhere valuesareimposedonquadraticboundaryinterpolantalongthenor- maldirection(reddashedline).Valuesforthebilinearinterpolants aresuppliedfromthegreencircledinteriorgridpoints........61 Figure2.9:tstudyfortheBesselmodeinacirculardomainwith CFLnumberof2.Usingquadraticboundaryinterpolantwithbilinear interiorinterpolant............................62 Figure3.1: S ( x )=1 j x j ; j x j < 1 ;S ( x )=0 ; j x j 1................71 Figure3.2:Thesupportofalinearparticleshape Sp ( x;y )in2D.........71 Figure3.3:Periodicextensionofshapefunctionforboundaryparticle......74 Figure3.4:Particleshapefunction S p ( x )ofaboundaryparticlewithghostcell [ x 0 ;x 1 ]=[ x; 0]............................82 Figure3.5:Particleshapefunctions S p ( x )(solid)and S img ( x )(dotted)ofabound- aryparticleanditsimageparticle,withghostcells[ x 1 ;x 0 ]=[ x; x ] and[ x 0 ;x 1 ]=[ x; 0]..........................82 vii Figure3.6:Potentialenergyincoldplasmaoscillation.Greenisthe1Dnumerical result,blueisthe2Dnumericalresult,andredisthepredictionof lineartheory.Weseetheplasmafrequencyisaccuratelyreproduced inoursimulations.............................82 Figure3.7:Growthofthemodewithmaximumgrowthrateinthetwostream instability,correspondingto k =3 : 06.Greenisthe1Dnumerical result,blueisthe2Dnumericalresult(measuredalongthecentral y =0slice),andredisthepredictionoflineartheory.Weseethe correctgrowthrateisreproducedinoursimulations..........83 Figure3.8: Weseeselectedparticlephasespaceplotsforthetwostreaminstability problem.Theleftcolumnisfromthe1Dsimulation,whiletherightcolumn isfromthe2Dsimulation,followingasliceofparticlesinitialized alongtheline y = L y = 2. .........................84 Figure3.9:Landaudampingofthelowestmode,correspondingto k =0 : 5.Green isthe1Dnumericalresult,blueisthe2Dnumericalresult(measured alongthecentral y =0slice),andredisthepredictionoflinearthe- ory.Weseethatthecorrectdecayrateisreproducedinoursimulations95 Figure3.10: In3.10a,weseethescalarpotentialat t =3 : 6thermal-iontransit times.In3.10b,weseethesimulationelectronandioncount,theredand bluecurvesrespectively,alongwiththeinjectionrate,theblackdashedline. 96 Figure3.11:ThestaggeredgridusedfortheBennettpinchproblem........97 Figure3.12: Theshows3.12aelectrondensity,3.12bandpotentialenergies,3.12c magneticatvariousCFLnumbers,and3.12dtherelativeerrorinthe azimuthalmagnetic B (normalizedtothemaximumvalueofthe magneticResultsarewithaCFLnumberof3,exceptasnotedin 3.12c.Positionunitsareintermsofeebeamradius R b . ......98 Figure3.13:ThestaggeredgridusedfortheMardahlbeamproblem........99 Figure3.14: Theshowsthedivergenceerrorintheelectricandthe beamdistributioncalculatedfromawaveequationpotential3.14a,3.14c andinthePoissonequationpotential3.14b,3.14d. ............100 viii Chapter1 Introduction Collisionlessplasmas-systemsofchargedparticlesinteractingthroughelectromagnetic -aremodelledbytheVlasov-Maxwellsystemofpartialtialequations(PDEs),which coupleMaxwell'sequations,describingtheevolutionoftheelectricandmagnetic E and B ,toVlasovequations,atypeofhyperbolicPDEdescribingtheevolutionofthephase-space distributionfunctions(DFs) f s ofthevariousspecies s ofchargedparticles.Particle-in-cell (PIC)methods[4,25,55],indevelopmentandusesincethe1960sandaprimarytoolin thecomputersimulationofplasmas,combineanEuleriandescriptionofthewitha LagrangiandescriptionoftheDFs;thatis,areevolvedonamesh,whileDFsare representedbymovingparticleswhosetrajectoriesarecharacteristicsofthecorresponding Vlasovequation.Thus,PICmethodsrequireamethodtocomputetheeldsonameshand amethodtocomputeparticletrajectories,aswellasinterpolationtoolstoprovidefortheir coupling.Thisworkfocusesonanewmethodforthecomputationofthealongwith associatedinterpolationtechniques. UndertheLorenzgaugecondition,Maxwell'sequationsreducetouncoupledwaveequa- tionsforthescalarandvectorpotentials,and A .Recently,anovelmethodforthesolution ofthewaveequationhasbeendeveloped[10,9,11],basedontheMethodofLinesTranspose (MOLT),dimensionalsplittingandant1Dintegralsolutionmethod,whichisun- conditionallystable(orA-stable)-thatis,itisnotsubjecttotheCourant-Friedrichs-Lewy (CFL)restrictionlimitingtheratioofthetemporalstepsizetothespatialstepsize,typical 1 ofwidelyusedexplicitmethods.Inthiswork,weapplythismethodtotheuncoupledwave equationsforand A tosolveMaxwell'sequationswithamethodcomparableincomputa- tionalcostandcomplexityofcodetoexplicitmethodssuchasthewell-knownYeescheme [56,53],butwithoutintroducingaCFLrestrictionbasedonthespeedoflightasinsuch explicitmethods.Inconjunctionwithanappropriatedescriptionofparticles,weseekto developaPICmethodthatretainsthesimplicityofexplicitmethods whileeliminatingthisCFLrestriction. Wedemonstratetheapplicationofourmethodtobothelectrostaticandelectromagnetic problems.Inthenon-relativistic,zero-magneticlimit,itistypicaltomaketheelec- trostaticapproximation, E = 2 = 0 ,simplifyingtheVlasov-Maxwellsystem totheVlasov-Poissonsystem.Correspondingly,inthisworkweconsiderthesamenon- relativistic,zero-magneticlimit,anddrop A andthecorrespondingwaveequations fromourmodel.Whenparticlevelocitiesaresmallcomparedtothespeedoflight,wear- guethatthismodel,whichwetermtheVlasov-Wavemodelduetothereplacementofthe Poissonequationwithawaveequation,willagreeapproximatelywiththeusualelectrostatic Vlasov-Poissonmodel.Wepresentnumericalresultsapplyingourmethodtoseveralstan- dardelectrostatictestproblemsinoneandtwodimensions,showingagreementwiththe predictionsoflineartheoryfortheelectrostaticmodel.Weapplyourmethodtoafully electromagneticbeampinchproblemintwodimensions.Itiswellknownthatanelectro- magneticPICmethodmustsatisfyadiscreteformofGauss'law( r E = 0 )toprevent seriousnumericalerrorsrelatedtotheviolationofchargeconservation[37].Inthiswork, weobtainsolutionsthatsatisfyexactlydiscreteformsofGauss'lawfortheelectric andthedivergence-freeconditionforthemagneticthroughastaggeredgridapproach, adaptedfromthewell-knownYeegrid[56],withaPoissonequationformulationforthescalar 2 potential.InadditiontoeliminatingtheCFLrestriction,thewavesolvermethodusedin thisworkthehandlingofcomplexboundarygeometryinaCartesiangridwithout usingastaircasingapproximation[9,54],andcanbeextendedtohigher-orderaccuracy[11], featuresthatwillbeincorporatedintoourPICmethodinfuturework. BuildingonpriorworkonunconditionallystableADI-FDTDschemesforMaxwell'sequa- tions[58,57,42,34],unconditionallystableADI-FDTDmethodsforMaxwell'sequations thatconservethedivergenceofthediscreteelectricandmagneticwereincorporated intoPICmethodsin[51],usingthemethodof[20]tohandlecomplexgeometries.Further, aPICmethodbasedonahigh-orderdiscontinuousGalerkinschemeforMaxwell'sequations wasdevelopedin[27,26],allowingforthehandlingofcomplexgeometriesthroughunstruc- turedmeshes,butwhichisstillsubjecttoCFLrestrictions.Futureworkwillcompareour methodtotheseapproachesinthecontextofPICmethods. InmostPICmethods,furtherstabilityrestrictionsalsoapply,themostrestrictiveoften beingtheneedtoresolvetheelectronplasmaperiod.Ourmethod,makinguseofexplicit algorithmsforadvancingparticles,issubjecttosuchrestriction.Inproblemswheretime scalesmuchlongerthantheelectronplasmaperiodareofprimaryinterestanddynamicson thescaleoftheelectronplasmaperiodcanbesafelyunderresolved-forinstance,incertain problemsinthestudyofiondynamics-itisdesirabletotakeamuchlargertimestepthan prescribedbytheplasmaperiodstabilityrestriction.Sinceyearsago[18,7,33,24],andwith arecentresurgence[13,38],ithasbeensoughttodevelopimplicitPICalgorithmsthatare notsubjecttothestabilityrestrictionbasedontheplasmaperiod(orthecyclotronperiod). Anultimatechallengeforourmethodwouldbesynthesiswithasuitableimplicitparticle integrationschemetoachieveapracticalfullyimplicitmethod,eliminatingthestability restrictionbasedontheplasmaperiodaswellastheCFLrestriction.However, 3 thatisbeyondthescopeofthepresentwork. 1.1Overview Wenowdescribethecontentofthisdissertation,andpointouttheoriginalcontributions giveninit.Thisdissertationbuildsontheimplicitwavesolverdevelopedin[10,9,11]. Thecontributionsofthepresentauthorto[9]includethedevelopmentofanembedded boundarymethodforNeumannboundaryconditionsonacomplexboundarygeometryand severalnumericalexamples,whicharereproducedinthisdissertation.Thedevelopmentof theembeddedboundarymethodforNeumannboundaryconditionsisimportant,asitis necessarytoimposeperfectelectricconducting(PEC)boundaryconditionswithourfor- mulationofMaxwell'sequations.Ananalysisin[28]suggeststhatsuchamethodwillbe unstablewithoutdissipation,sothepresentauthordevelopedamethodfor dissipationincenteredschemesthatmaintainsunconditionalstablity.Furthercontributions ofthepresentauthorconsistofthedevelopementofaPICmethodforthesimulationof plasmasbasedontheimplicitwavesolver,whichisgiveninthisdissertation.Important aspectsofthisPICmethodareafastconvolutionalgorithmforparticlesourcetermsanda staggered-gridapproachtoenforcingdivergenceconditions.Wenowdescribethestructure ofthepresentworkwithrespecttothepreviousworkontheimplicitwavesolverandthe originalcontributionsofthepresentauthor. InSection1.2,wedescribethemathematicalmodelsusedinthePICmethod,andin Section1.3,wedescribethenondimensionalizationsusedinthiswork.InSection2.1,we giveaheuristicexplanationoftheCFLrestrictionintypicalexplicitwavesolversalongwith amethodtomitigatethisrestriction,andinSection2.2,wedescribethederivationoftwo 4 second-orderaccuratesemi-discreteschemes,describingworkfrom[11]and[9].InSection 2.3,wegiveasemi-discreteVonNeumannanalysisoftheseschemes,whichisanoriginal contibution.InSection2.4,wedescribethefastnumericalmethodusedforanintegral solution,describingworkfrom[9]. Theremainderofthedissertationisoriginalworkofthepresentauthor.InSection2.5, wegiveaproofthatthefullydiscreteschemeisunconditionallystable,inthesenseofVon Neumannanalysis.InSection2.6,wedescribeamethodforintroducingdissipation intocenteredschemesofarbitraryorderaccuracywhilemaintainingunconditionalstability. InSection2.7,wedescribeafourth-orderaccurateNewmarkmethodforthewaveequation thatpossessesdissipation.InSection2.8,wedescribeanembeddedboundarymethodfor thewaveequationwithNeumannboundaryconditionsonacomplexboundarygeometry.In Table2.1,wegiveatablethatsummarizestheimplicitwavesolver(includingtheprevious workdescribedabove).InSection2.9,wepresentseveralnumericalexamplesoftheimplicit waveequationsolver. ThePICmethodisdescribedinChapter3.Thechargeandcurrentweightingschemes usedaregiveninSection3.1,andamethodforreducingdivergenceerrortomachineprecision isgiveninSection3.2.TheparticleequationsofmotionaredescribedinSection3.3.A fastalgorithmforevaluatingconvolutionsofparticlesourcetermsisgiveninSection3.4. SomeissueswithimposingboundaryconditionswithparticlesaredescribedinSection3.5. InSection3.6,wepresentseveralnumericalexamplesusingourPICmethod.Finally,we concludethedissertationwithabriefdiscussioninChapter4. 5 1.2MathematicalModels Theself-consistentevolutionofacollisionless,singlespeciesplasmaisdescribedbythe Vlasov-Maxwellsystem,givenintheSIsystemofunitsby @f @t + v r x f + q m ( E + v B ) r v f =0(1.1) 0 0 @ E @t = r B 0 J @ B @t = E (1.2) r E = ˆ=" 0 r B =0(1.3) ˆ ( x ;t )= q Z v f ( x ; v ;t ) d v ; J = q Z v v f ( x ; v ;t ) d v (1.4) where f ( x ; v ;t )isthephase-spacedistributionfunction, q isthechargeand m isthemassofa particle, E ( x ;t )istheelectric B ( x ;t )isthemagnetic ˆ ( x ;t )isthechargedensity (and ˆ B representsastaticuniformbackgroundchargedistribution), J ( x ;t )isthecurrent density, " 0 istheelectricpermittivityofthevacuum,and 0 isthemagneticpermability ofthevacuum.(Boldfacevariablesaretostandforvectorquantities,whilenonboldface variablesaretostandforscalarquantities.) Inthelimit,inanappropriatesense,of j v j =c ! 0(where c =1 = p " 0 0 isthespeedof lightinvacuum)andintheabsenceofmagnetics,theVlasov-Maxwellsystemreduces 6 totheVlasov-Poissonsystem: @f @t + v r x f q m r r v f =0(1.5) = ˆ=" 0 ;ˆ ( x ;t )= q Z v f ( x ; v ;t ) d v (1.6) TheVlasov-MaxwellandVlasov-Poissonmodelshavebeenwidelystudied,andmany numericalmethodshavebeendevelopedtosolvethem,suchaselectrostaticandelectromag- neticPICmethods[4,25,17,39,23,27,51],aswellasEulerianmethods[22,45,1,16,15] andsemi-Lagrangianmethods[14,52,47,48].Insteadofsolvingeitherofthesesystemsdi- rectly,weseektodevelopasemi-implicitapproachforthebasedonavectorpotential formulation.UndertheLorenzgaugecondition, 1 c 2 @ @t + r A =0 ; Maxwell'sequationsreducetouncoupledwaveequations: 1 c 2 @ 2 @t 2 = ˆ=" 0 (1.7) 1 c 2 @ 2 A @t 2 A = 0 J (1.8) with E = @ A @t and B = r A [31].Asimilarpotential-basedapproachwas takenin[40],whichusedtheCoulombgaugecondition( r A =0),resultinginmore acomplicatedsetofcoupledequations,whichtheysolvedusingexplicit methods.WechoosetoworkinseadintheLorenzguage,astheresultinguncoupledwave 7 equationsaresimplertohandlenumerically.Wenotethatperfectelectricconducting(PEC) boundaryconditionsmaybeimposedonthepotentialsbysetting=0, A n = 0 and ( r ( A n )) n =0ontheboundary,where n denotesthenormaltotheboundary. Intheelectromagneticcase,wewillsolvethesewaveequationscoupledintoaparticle- in-cellmethodtoconstituteasolutionoftheusualVlasov-Maxwellsystem.Inthenon- relativistic,zero-magneticlimit,wedrop A andthecorrespondingwaveequationsand consideraquasi-electrostaticVlasov-Wavemodel,whichwillagreewiththeelectrostatic Vlasov-Poissonmodelwhen j v j =c<< 1.Thisconditionfrequentlycanbeinterpretedas ! p L<T +2 R=cT ,wehave u ( x;t )= O (( ct ) 1 ) d =2(1.30) u ( x;t )=0 d =3(1.31) Asageneralizationofthis,foranytlysmooth f ( ;t )supportedin B (0 ;R )forall t> 0with f ( ;t )= f T ( )forall t>T ,itisagaineasilyarguedthatfor x 2 B (0 ;R )and 11 t>T +2 R=cT ,wehave u ( x;t )= u 2 P ( x )+ O (( ct ) 1 ) d =2(1.32) u ( x;t )= u 3 P ( x ) d =3(1.33) where u d P ( x )istheclassicalintegralsolutionofthePoissonequation u d P = f T in dimension d . TheconvergenceofsolutionstotheVlasov-MaxwellsystemtothoseoftheVlasov-Poisson systemhasbeenrigorouslyconsideredinworkssuchas[2,50].Itmaybepossibletoapply similartechniquestorigorouslystudytheconvergenceofsolutionsofourVlasov-Wavemodel tothoseoftheVlasov-Poissonsystem,butthisisoutsideofthescopeofthepresentwork. 1.3.2ElectromagneticCase HerewegivethenormalizationusedfortheelectromagnetictestprobleminSection3.6.2. Considerthefollowingchangeofvariables: ~ f = Ff; ~ t = Tt; ~ x = L x ; ~ v = V v ; (1.34) ~ ˆ = qNˆ; ~ J z = qVNJ z ; (1.35) ~ = 0 ; ~ A = A 0 A z (1.36) appliedtothesystem @ ~ f @ ~ t + ~ v r ~ x ~ f + q m ~ x ~ ~ v ( r ~ x (0 ; 0 ; ~ A z )) r ~ v ~ f =0(1.37) 12 1 c 2 @ 2 ~ @ ~ t 2 ~ =~ ˆ=" 0 ; ~ ˆ ( ~ x ; ~ t )= q Z ~ v ~ f ( ~ x ; ~ v ; ~ t ) d ~ v (1.38) 1 c 2 @ 2 ~ A z @ ~ t 2 ~ A z = 0 ~ J z ; ~ J z ( ~ x ; ~ t )= q Z ~ v ~ v z ~ f ( ~ x ; ~ v ; ~ t ) d ~ v (1.39) Assumingthescalings, V = L T ;T = r m" 0 Nq 2 = ! 1 p ; (1.40) 0 = qNL 2 " R " 0 ;A 0 = 0 VqNL 2 R ;" R R = V 2 =c 2 ; (1.41) F = N 2 d d ( " 0 m ) d= 2 q d L d (1.42) (1.43) where " R and R aredimensionlessparametersintroducedtoenforcetheLorenzgauge condition,weobtain: @f @t + v r x f + " R ( + v ( r (0 ; 0 ;A z )) r v f =0(1.44) 2 @ 2 @t 2 = ˆ=" R ;ˆ ( x ;t )= Z v f ( x ; v ;t ) d v (1.45) 2 @ 2 A z @t 2 A z = R J z ;J z ( x ;t )= Z v v z f ( x ; v ;t ) d v (1.46) where = L cT = V c ,asintheelectrostaticcase.FortheBennettpinchprobleminSection 3.6.2,wechoose " R =1. 13 Chapter2 Descriptionofthewaveequation solver 2.1OriginandMitigationoftheCFLRestriction ConsidertheCauchyproblemforthewaveequation: u tt = c 2 u (2.1) u ( x; 0)= g ( x )(2.2) u t ( x; 0)= h ( x ) : (2.3) ToseetheoriginoftheCFLrestrictioninatypicalexplicitmethod, considerthesemi-discreteschemeobtainedbysubstitutingthecentered,second-order discretizationof u tt intothewaveequation[19].Weobtain u ( x;t + t ) 2 u ( x;t )+ u ( x;t t ) t 2 = c 2 2 u ( x;t ) : UponFouriertransforminginthespatialvariable,weobtain 14 ^ u ( ˘;t + t ) 2^ u ( ˘;t )+^ u ( ˘;t t )= c 2 t 2 j ˘ j 2 ^ u ( ˘;t )(2.4) where ˘ isthespatialfrequency.Theproblemisthatahigh-frequencyperturbationof u -suchasthatintroducedbytruncationerror-willbebythesymbol c 2 t 2 j ˘ j 2 , resultingininstability.Inafullydiscretemethod,thespatialfrequenciesareboundedby 2 = x (theNyquistfrequency),meaningthat c t= x mustbesutlysmalltoprevent ofhighfrequencies.Intermsofsemi-discreteVonNeumannanalysis,consider substitutinginto2.4asemi-discretesolutionoftheform^ u ( ˘;n t )= n forsome 2 C . Weobtainaquadraticequation 2 (2 z 2 ) +1=0andthestabilitycondition j j 1 foranyroot ofthisequation,where z = c t j ˘ j .Itiseasilyvfromthequadratic formulathattheroots i satisfy j i j 1for i =1 ; 2ifandonlyif j z j 2,thatisifandonly if c t 2 = j ˘ j .Forthe1Dcase,insertingtheNyquistfrequency2 = x intotheupperbound givestheusualCFLstabilityrestriction c t x . ToavoidtheCFLrestrictionwhilemaintaininganexplicitupdateequation,wemodify thesymboloftheLaplaciantobounditandpreventtheamplofhighfrequencies. Considertheone-dimensionalcase.Welookforboundedrationalapproximationsof ˘ 2 near ˘ =0.Oneoptionis ˘ 2 ˇ 2 ˘ 2 ˘ 2 + 2 ˇ ˘ 2 (1 ( ˘ ) 2 + ) Wecanwritethismosemi-discreteschemeinFourierspaceas 15 ^ u ( ˘;t + t ) 2^ u ( ˘;t )+^ u ( ˘;t t )= c 2 t 2 2 ˘ 2 ˘ 2 + 2 ^ u ( ˘;t ) Choosing = =c t forsomeproperlychosen ,wecanseethatthemosymbol isuniformlyboundedby 2 forall t ,sothatchoosing 2willensurestability,inde- pendentlyof t ,anditturnsoutthat =2givesoptimalaccuracyforstableschemesof thisfamily[9,11].MultiplicationinFourierspaceisequivalenttoconvolutioninphysical space,sotransformingbackintophysicalspacegivesamethodbasedonconvolutionwith theinverseFouriertransformofthemosymbol: u ( x;t + t ) 2 u ( x;t )+ u ( x;t t )= 2 D [ u ]( x;t ) where D [ f ]( x )= R ( ( x x 0 ) 2 e j x x 0 j ) f ( x 0 ) dx 0 = f ( x ) 2 R e j x x 0 j f ( x 0 ) dx 0 . Asageneralcomment,thismoofthesymboloftheLaplacianintroduceserror athighspatialfrequencies,sothatthistechniquemaynotbesuitableforproblemsinwhich extremelyaccuratepropagationofwavesathighspatialfrequencies(highwavenumber)is essential.However,inmanyproblemsofplasmaphysics,thephysicsisdominatedby atlowspatialfrequencies,andintheseproblemsthistechniquecanusefulapplication. Adetaileddiscussionofnumericalmethodsforthewaveequationderivedfromtheper- spectiveofFouriermultipliers,includinghigher-ordermethodsandtheextensiontomultiple dimensions,canbefoundin[11].Inthefollowingsection,wegivethederivationoftwouseful second-ordernumericalmethodsthroughalternativemeans,theMethodofLinesTranspose, 16 whichalsofacilitatestheextensionofthemethodstoboundeddomainsandtomultiple dimensions. 2.2DerivationofSecond-OrderSemi-DiscreteSchemes WiththeMethodofLinesTranspose IntheMethodofLinesTransposeforthesolutionoftime-dependentPDEs,alsoknown asRothe'smethod[49],discretizationsoftimederivativesaresubstituted intothePDE,resultinginaboundaryvalueproblem(BVP)tobesolvedateachtime step.Inrecentyears,theMOLThasbeenappliedasanumericalmethodtosolvevarious time-dependentPDEs,withinitialfocusonparabolicequations[12,29,30,8,36].The work[10]extendedthisapproachtothesecondorderwaveequation,andisthebasisofthe presentwork.Furtherworkextendedthisapproachtohigherdimensionsthroughdimen- sionalsplitting[9]andtohigher-orderthroughtheFouriermultiplierapproachmentioned intheprevioussection[11].Basedontheworkin[10,9,11],wegivehereanoverviewof thederivationoftwousefulsecond-orderschemesforthewaveequation,theirextensionto multipledimensionsthroughdimensionalsplitting,andafastnumericalalgorithmforthe 1Dproblem.InSection2.3,wegivesemi-discreteVonNeumannanalysesanddispersion relationsforthesesemi-discreteschemes,andinSection2.5wegiveanewproofthatthe fullydiscreteschemesareunconditionallystable,inthesenseofVonNeumannanalysis. 17 2.2.1DispersiveSemi-DiscreteScheme Asintheprevioussection,wesubstitutethefollowingsecond-ordercentereddiscretization intothewaveequation 1 c 2 u tt u = f ( x;t ): u n tt = u n +1 2 u n + u n 1 t 2 t 2 12 u tttt ( x; ) : InsteadofprocedingbytheFouriermultiplierapproachasmentionedpreviously,we applyanaveragingtechniquewithsimilarresults.WeaveragetheLaplaciantermintime u n = 1 4 u n +1 +2 u n + u n 1 + O t 2 ) andsubstituteintothewaveequation.ng =2 = ( c t )andmanipulatinggivesthe semi-discretescheme 1 2 +1 u n +1 +2 u n + u n 1 =4 u n + 4 2 f ( x;t n )+ O t 4 ) : Wecallthisthedispersivesemi-discretescheme,sincealltermsinthesemi-discrete dispersionrelationarereal-valued(see2.3). 2.2.2eSemi-DiscreteScheme Wesubstitutethefollowingbackwardformula(BDF)discretization: u n +1 tt = 2 u n +1 5 u n +4 u n 1 u n 2 t 2 t 2 12 u tttt ( x; ) intothewaveequation 1 c 2 u tt u = f ( x;t ). 18 Rearranging, = p 2 = ( c t )anddividingby 2 givesthesemi-discretescheme 1 2 +1 u n +1 = 1 2 5 u n 4 u n 1 + u n 2 + 1 2 f ( x;t n +1 )+ O t 4 ) : Wecallthistheesemi-discretescheme,duetothepresenceofimaginary-valued termsinthesemi-discretedispersionrelation(see2.3).The(slight)dissipationofthismethod provedusefulininitialPICsimulationsaimedatcapturingthesteady-statepotentialdueto asingleparticle,andsoweusethisscheme(ratherthanthedispersivescheme)inallPIC simulationsinthiswork.Ontheotherhand,thedispersivescheme(andthehigher-order versionsin[11])havesmallertruncationerrorthantheescheme,anddonothavean implicitsourceterm(attimelevel n +1)thatconfoundtheuseofhigher-ordertimestepping schemesforparticles.Futureworkwillinvestigatetheuseofthesetime-centeredschemes, makinguseofdissipation(asinSection2.6)whenappropriate,andcorresponding time-steppingschemesforparticles. 2.3Semi-DiscreteVonNeumannAnalysisoftheSec- ondOrderSchemes Inthissection,weprovidesemi-discreteVonNeumannstabilityanalysesanddispersion relationsforthesemi-discreteschemesderivedin2.Thesebuildonsimilaranalysesfor related,buterent,schemesgivenin[10]. 19 2.3.1DispersiveScheme Substituteansatz u n = e i ( k x ~ !n t ) into 1 2 +1 u n +1 +2 u n + u n 1 =4 u n Obtainapolynomial 2 +2 j k j 2 2 j k j 2 + 2 +1=0 where = e i ~ ! t . Wecansolvetoobtain 1 = 2 = i j k j i j k j ,whichgives j 1 = 2 j =1,meaningthisschemeis non-dissipative. Notingthatcos(~ ! t )= 1 2 ( 1 + 2 ),and z = j k j = j k j c t= 2,weobtain ~ ! = 2 t arccos s 2 j k j 2 + 2 ! = 2 t arccos r 1 1+ z 2 ! ˇ 2 t z z 3 = 3+ O ( z 5 ) ˇj k j c 1 1 12 ( j k j c t ) 2 + O (( j k j c t ) 4 ) Sothephaseerroris j ~ ! j k j c 1 j = O (( j k j c t ) 2 ).Moreover,~ ! isreal,owingtothenon- dissipativenatureofthescheme. 2.3.2eScheme Substitutingtheansatz u n = e i ( k x ~ !n t ) into 1 2 +1 u n +1 = 1 2 5 u n 4 u n 1 + u n 2 ; 20 weobtainapolynomial 3 4 2 +5 2(1+ z 2 )=0 where = e i ~ ! t , z = j k j .Thethreerootstellusaboutpossiblemodes u n = e ik x n . Anecessaryconditionforstabilityis 1forallroots. Theroot, 1 = 4 3 + 1 3 27 z 2 +3 p 3 p 27 z 4 +2 z 2 +1 1 = 3 + 1 3 27 z 2 +3 p 3 p 27 z 4 +2 z 2 +1 1 = 3 ; correspondstoaspuriousnonpropagatingmodeoftheform u n = e ik x n 1 .Since 1 2 forall z = j k j c t= p 2 0,themoderapidlydecaysandposesnothreattostability. Theothertworootsareapairofcomplexconjugates: 2 = 3 = 1 6 (1 i p 3) 27 z 2 +3 p 3 p 27 z 4 +2 z 2 +1 1 = 3 1 i p 3 6(27 z 2 +3 p 3 p 27 z 4 +2 z 2 +1) 1 = 3 +4 = 3 ˇ 1 i p 2 z z 2 i 5 p 2 4 z 3 +4 z 4 i 231 p 2 32 z 5 + O ( z 6 )as z ! 0 : Wecanshowthat j 2 = 3 j 2 = 16 9 + 4 W 4 16 W 3 4 W 2 16 W +4 36 W 2 1forall W 1 21 where W = 27 z 2 +3 p 3 p 27 z 4 +2 z 2 +1 1 = 3 1forall z 0. Since = e i ~ ! t ,wehave ~ ! = 1 i t log( ) ˇ 1 i t log(1+ i p 2 z z 2 i 5 p 2 4 z 3 +4 z 4 + i 231 p 2 32 z 5 + O ( z 6 )) ˇj k j c 1 11( j k j c t ) 2 24 i ( j k j c t ) 3 2 15( j k j c t ) 4 16 + O (( j k j c t ) 5 ) Sothephaseerroris j ~ ! j k j c 1 j = O (( j k j c t ) 2 ). Thepresenceoftheimaginarythirdtermintheexpansionshowsthat u n ˘ e in ( i ( j k j c t ) 3 t ) = e n ( j k j c ) 3 t 4 = 2 ,causingthemodetodecay.Thisiswhywetermthisschemee(or dissipative). 2.4SolutionoftheMoHelmholtzEquation Asseenintheprevioussection,boththedispersiveandesemi-discreteschemes includeanellipticBVPtobesolvedateachtimestep.TheresultingPDEissometimes calledthemoHelmholtzequation[30].Incontrasttotheusualfrequency-domain Helmholtzequation u + ! 2 c 2 u = f ,themodHelmholzequationhasanonoscillatory Green'sfunction.Theoscillationinthesolutionofthewaveequationissupportedbythe presenceofmultipletimelevelsinthesemi-discreteequations.Oursolutionstrategyisto usethewell-knowntechniqueofdimensionalsplitting[43]toreduceproblemsinmultiple dimensionstoproblemsinonedimension,towhichweapplyafastintegralsolutionmethod. Thedimensionally-splitintegralsolutionnaturallyleadstounconditionallystablenumerical 22 schemeswithcomputationalcostandcodingcomplexitycomparabletoexplicitschemes. 2.4.1DimensionalSplitting Forsmooth u ( x;y ),wehave 1 2 @ xx + @ yy +1 u = 1 2 @ xx +1 1 2 @ yy +1 u 1 4 @ xxyy u So,toapproximatelysolve 1 2 @ xx + @ yy +1 u = f wewillinsteadsolve 1 2 @ xx +1 1 2 @ yy +1 u = f: wherewehaveintroducedthesplittingerror 1 4 @ xxyy u = O ( c 4 t 4 ) : Tosolve 1 2 @ xx +1 1 2 @ yy +1 u = f we w = 1 2 @ yy +1 u andsolvethefollowingone-dimensionalBVPs\linebyline": 1 2 @ xx +1 w ( x; )= f ( x; ) 1 2 @ yy +1 u ( ;y )= w ( ;y ) whereappropriateboundaryconditionsaresupplied(see[10,11]).Since w = u + O ( c 2 t 2 ), forsecondorderaccuracyittousethesameboundaryconditionsfor w as u .Future 23 workwillinvestigatehigher-orderaccurateboundaryconditions.Tofacilitatethe\line-by- line"solution,wediscretizethedomainwithauniformCartesiangrid. 2.4.2IntegralSolutionMethodin1D ConsiderthemoHelmholtzequationinonedimension, 1 2 d 2 dx 2 +1 u = f ( x )for x 2 =( a;b ),withappropriateboundaryconditionsimposedat x = a and x = b .Wecan write u ( x )= u P ( x )+ u H ( x ),where u P ( x )= 2 Z b a f ( x 0 ) e j x x 0 j dx 0 u H ( x )= Ae ( x a ) + Be ( b x ) aretheparticularandhomogeneoussolutionsandwhere A and B dependonthebound- aryconditionsimposed,aswellasthevaluesof u P ( a )and u P ( b )[10].Ourfastintegral solverconsistsofafastconvolutionalgorithmfortheevaluationoftheparticularsolution u P ,alongwithappropriatealgorithmsforevaluatingthehomogeneoussolution u H ,which canbeviewedasboundarycorrectionterms.Formanycommonboundaryconditions,the cots A and B fortheboundarycorrectiontermscanbefoundbyapplyingthegiven boundaryconditionsandsolvinga2 2systembyhandfor A and B .Forinstance,inthecase ofhomogeneousDirichletboundaryconditions u ( a )= u ( b )=0and = e ( b a ) , we A = u P ( b ) u P ( a ) 1 2 B = u P ( a ) u P ( b ) 1 2 : 24 Inthecaseofperiodicboundaryconditions,we A = u P ( b ) 1 B = u P ( a ) 1 : Withsomefurtherconsideration,otherboundaryconditionscanbederived,including w(absorbing)boundaryconditions(fortheunderlyingwaveequation,basedonone- waywaveequations).Forfurtherdetails,see[10,9,11]. 2.4.3FastNumericalEvaluationofthe1DConvolutionOperator Theconvolutionoperatorgivingtheparticularsolutioncanbedecomposedas I [ f ]( x )= 2 Z 1 f ( x 0 ) e j x x 0 j dx 0 = 2 Z x f ( x 0 ) e j x x 0 j dx 0 + 2 Z 1 x f ( x 0 ) e j x x 0 j dx 0 =: I L [ f ]( x )+ I R [ f ]( x ) foragivenfunction f .Meanwhile,wehavetherecursionrelations I L [ f ]( x )= e x I L [ f ]( x x )+ 2 Z x x x f ( x 0 ) e ( x x 0 ) dx 0 I R [ f ]( x )= e x I R [ f ]( x + x )+ 2 Z x x x f ( x 0 ) e ( x 0 x ) dx 0 : Basedontheseobservations,weoutlinethefastalgorithmforthenumericalevaluation ofthisconvolutionoperatoronauniformCartesiangriddevelopedin[10,9,11].Consider theconvolutionoperatorappliedtoafunction f supportedontheinterval( a;b ),withthe 25 convolutionalsotobeevaluatedin( a;b ).Theintervalisdiscretizedinto N equalsubintervals oflength x =( b a ) =N ,withendpoints x 1 = a , x j +1 = x j + x for j =1 ;:::;N .We denote I j = I [ f ]( x j ), I L j = I L [ f ]( x j )and I R j = I R [ f ]( x j ),asabove.Further,we thelocalintegrals J L j = 2 Z x j x j 1 f ( x 0 ) e ( x j x 0 ) dx 0 j =2 ;:::;N +1 J R j = 2 Z x j +1 x j f ( x 0 ) e ( x 0 x j ) dx 0 j =1 ;:::;N: Supposewehavealreadycomputedthese J L j and J R j .Setting I L 1 = I R N +1 =0,we canthenperformarecursiveevaluationoftheconvolutionintegralbycomputing I L j +1 = e x I L j + J L j +1 and I R N +1 j = e x I R N +1 j +1 + J R N +1 j for j =1 ;:::N .Wethensum I j = I L j + I R j foreach j =1 ;::N +1. Themethodforevaluatingthelocalintegrals J L j and J R j dependsonthenatureofthe integrand.Fortheparticleconvolutionintegral,wheretheintegrandisthesumofparticle shapefunctions,wecananalyticallyevaluatethelocalintegrals.ThisisdescribedinSection 3.4.Forgeneralintegrands,andspforthetermsinvolving u ontherighthandsides ofthesemi-discreteschemes,wenumericallyevaluatethelocalintegralswithaquadrature rule.Itisimportanttonotethatanygivenquadraturerulemayormaynotdeliveran accurateandstableoverallschemeforthewaveequation.Toachieveaccuracyandstability, weuseaquadraturerulefoundbyanalyticallyintegratingagainstaLagrangepolynomial interpolant.Foraquadraticinterpolant,leadingtosecond-orderaccuratequadraturein 26 space,weobtainthefollowingapproximations J L j ˇ Pf ( x j )+ Qf ( x j 1 )+ R ( f ( x j +1 ) 2 f ( x j )+ f ( x j 1 )) J R j ˇ Pf ( x j )+ Qf ( x j +1 )+ R ( f ( x j +1 ) 2 f ( x j )+ f ( x j 1 )) where, = x and d = e , P =1 1 d Q = d + 1 d R = 1 d 2 1+ d 2 Higher-orderspatialaccuracycanbeobtainedbyusinghigher-orderaccuracyquadrature rules,withtheoutcomeofunconditionalstabilitylimitingthechoiceofquadraturerules. Furtherdetailscanbefoundin[10,9,11].InSection2.5,VonNeumannanalysesarecarried outforthefullydiscreteanddispersiveschemesinoneandmultipledimensions,and itisshownthattheyarebothunconditionallystable. 2.5FullyDiscreteVonNeumannAnalysis Inthissection,weprovidefullydiscreteVonNeumananalysesforthetwofullydiscrete schemesderivedin2,andshowthattheyareunconditionallystable,inthesenseofVon Neumannanalysis.Combiningthequadraturerulesandexponentialrecursion,andignoring 27 boundaries,wecanwrite I [ f ]( x j )= 2 Z 1 f ( x 0 ) e j x j x 0 j dx 0 ˇ ˇ I h [ f j ]= Pf j + 1 2 Q f j +1 + f j 1 + R f j +1 2 f j + f j 1 + + 1 2 1 X k =1 e k P f j + k + f j k + Q f j + k +1 + f j k 1 + + R f j + k +1 2 f j + k + f j + k 1 + f j k +1 2 f j k + f j k 1 with = h = x and P , Q ,and R asinSection2. Usingthediscreteconvolutionoperator I h ,andignoringsourcesandboundaries,wecan writetheeversionofthefullydiscreteschemeas u n +1 j = 1 2 I h [5 u n j 4 u n 1 j + u n 2 j ] andthedispersiveschemeas u n +1 j +2 u n j + u n 1 j =4 I h [ u n j ] I h;x and I h;y asthediscreteconvolutionoperatorsactinginthe x -and y - directions,andagainignoringsourcesandboundaries,wecanwritetheeversionof the2Dfullydiscreteschemeas u n +1 jk = 1 2 I h;x [ I h;y [5 u n jk 4 u n 1 jk + u n 2 jk ]] 28 andthe2Dfullydispersiveschemeas u n +1 jk +2 u n jk + u n 1 jk =4 I h;x [ I h;y [ u n jk ]] Wewillrefertotheseasfree-spaceschemes.Wenowstateastabilitytheorembasedon theVonNeumannanalysisoftheschemes. Theorem. Thefullydiscretedispersiveandfree-spaceschemesinonedimen- sionareunconditionallystable.Inhigherdimensions,thecorrespondingdimensionally-split schemesarealsounconditionallystable. Toprovethestabilitytheorem,weconsidersomepropertiesofthediscretecovolution operator. Lemma. theationfactor A ( ~ k; )= e i ~ kj x I h [ e i ~ kj x ] .Then A ( ~ k; ) iswell- d(independentof j ),andthefollowinghold. If > 0 and 0 < j ~ k x j ˇ ,then 0 0 ,then A (0 ; )=1 . Forany 0 < j ~ k x j ˇ , lim ! 0 + A ( ~ k; )=0 . Forany 0 < j ~ k x j ˇ , lim !1 A ( ~ k; )=1 . ProofoftheLemma. Wecalculate: 29 A ( ~ k; )= e i ~ kj x I h [ e i ~ kj x ] = P 1+ 1 X k =1 e k cos( k ~ k x ) ! + + Q e 1 X k =1 e k cos( k ~ k x ) ! + +2 R (cos( ~ k x ) 1) 1+ 1 X k =1 e k cos( k ~ k x ) ! =1+ T; T = d 2 + d (cos( ~ k x ) 1) 2 2 1 d 2 ( d cos( ~ k x ) 1)(cos( ~ k x ) 1) d 2 2cos( ~ k x ) d +1 with d = e . If ~ k x =0and > 0,then T =0.If0 < j ~ k x j ˇ ,somecalculusshows 1 0,andlim ! 0 + T = 1,andlim !1 T =0. ProofoftheTheorem. Weprovethateachone-dimensionalfree-spacefullydiscrete schemeisunconditionallystable,thendescribetheextensionoftheresulttomultipledi- mensionsforthedimensionally-splitschemes. DispersiveScheme Inthissection,weprovethattheone-dimensionalfree-spacefullydiscretedispersivescheme isunconditionallystableandnon-dissipative.Considerapplyingthefullydiscretedispersive schemeto u n j = e i ( ~ kj x ~ !n t ) .Weobtainapolynomial 2 +(2 4 z ) +1=0 30 where = e i ~ ! t and z = A ( ~ k; ).Notethatthelemmaimplies0 z 1forall ~ k and x; t> 0.Wecansolvetoobtaintheroots 1 = 2 =(2 z 1) i p (1 z ) z .Wecanshow j 1 = 2 j =1forall0 z 1,sothefullydiscretedispersiveschemeisunconditionallystable, andnon-dissipative.Thisprovesthetheoreminthecaseoftheonedimensionaldispersive scheme. eScheme Considerapplyingthefullydiscreteeschemeto u n j = e i ( ~ kj x ~ !n t ) .Weobtaina polynomial 3 4 2 +5 z =0 where = e i ~ ! t and z =2 =A ( ~ k; ).Notethatthelemmaimplies z 2forall ~ k and x; t> 0.Theconditionforstabilityis j j 1forallrootsofthepolynomial,whichwe verifybelow. Therootcorrespondstoaspuriousmode: 1 = 4 3 + 3 p 3 p 3 p 27 z 2 104 z +100+27 z 52 3 3 p 2 + + 3 p 2 3 3 p 3 p 3 p 27 z 2 104 z +100+27 z 52 Wecanshowthat 1 = W 2 +4 W +1 3 W 2forall W 1 where W = 3 q 3 p 3 p 27 z 2 104 z +100+27 z 52 3 p 2 1for z 2. 31 Theotherrootsarepairofcomplexconjugates: 2 = 3 = 4 3 (1 i p 3) 3 p 3 p 3 p 27 z 2 104 z +100+27 z 52 6 3 p 2 (1 i p 3) 3 p 2 6 3 p 3 p 3 p 27 z 2 104 z +100+27 z 52 Wecanshowthat j 2 = 3 j 2 = 4 W 4 16 W 3 +60 W 2 16 W +4 36 W 2 1forall W 1 where W = 3 q 3 p 3 p 27 z 2 104 z +100+27 z 52 3 p 2 1for z 2. Thisprovesthetheoreminthecaseoftheonedimensionalescheme. ExtensiontoHigherDimensions Whenapplyingthedimensionally-splittwo-dimensionalschemesto u n jk = e i ( ~ k x j x + ~ k y k y ~ !n t ) , weobtainthesameVonNeumannpolynomialsasinthe1Dcase,andbasicallythesamesta- bilityanalysiscanberepeated.Thiscanbeeasilyextendedtodimensionally-splitschemes inanydimension. ThisVonNeumannanalysisdoesnottakeintoconsiderationtheofboundary conditions,andinprinciple,certainnumericalboundaryconditionscouldresultininstability. Inthetestproblemspresentedin[10,9,11]aswellasthepresentwork,thestabilityof themethodseemsrobustunderavarietyofnumericalboundaryconditions.Astability analysisforsome1Dfullydiscreteschemes(slightlytfromthosepresentedhere) 32 withnumericalDirichletandNeumannboundaryconditionswascarriedoutin[10],showing unconditionalstabilityinthoseschemes.Asimilarstabilityanalysiscouldalsobecarried outfortheschemesconsideredinthisworktostudythestabilityoftheseschemesunderthe inclusionofnumericalboundaryconditions. 2.6DissipationinCenteredSchemes Asdiscussedin2.8,itisnecessarytoincludesomedissipationinthenumerical schemetomaintainstabilitywithembeddedboundarymethodsforNeumannboundary conditions,necessaryfortheimplemenationofPECboundaryconditionswithcomplexge- ometries.Theeschemeisdissipative,buthasanimplicitsourceterm(attimelevel n +1),presentingayincouplingwithparticlesinthefullyelectromagneticcase. Thisseemstorestrictthemethodbasedontheesolvertouseaccurate, explicitforwardEulertimesteppingorsecond-orderaccurate,implicitCrank-Nicholsontime steppingforparticles,thesecondoptionrequiringiterationoftheandparticleupdates. Hence,weseektoschemesthatmaintainexplicitsourcetermsbutalsoincludesome dissipation.Thesecond-orderdispersiveschemegivenabove,andthehigher-orderaccurate schemesgivenin[11],haveexplicitsourceterms(attimelevel n ),butarenon-dissipative. Inthissection,wegiveamethodforaddingtunableialdissipationtermsintothese second-andfourth-orderaccuratecenteredschemes. 2.6.1Dissipationin1D Wegivemoformsofthecenteredversionsoftheimplicitwavesolverwithtunable dissipationthatretainthepropertyofunconditionalstability.Werestrictour 33 attentiontothe1Dschemeshere,butthesameproofsshouldapplyinthehigher-dimensional case.Welet denoteasmalldissipationparameter,and D [ u ]= u L 1 [ u ]= u ( x ) 2 R 1 e j x x 0 j u ( x 0 ) dx 0 bedasusual. Motheschemesgivenin[11],wehavethesecondorderschemewithdissipation, u n +1 2 u n + u n 1 = 2 D [ u n ]+ D 2 [ u n 1 ] ; (2.5) andthefourthorderschemewithdissipation, u n +1 2 u n + u n 1 = 2 D [ u n ] ( 2 + 4 12 ) D 2 [ u n ]+ D 3 [ u n 1 ] : (2.6) Wenowprovetheunconditionalstabilityoftheseschemesforprescribedvaluesof .As in[11],wepasstothehigh-frequencylimit. SecondOrderScheme WeobtaintheVonNeumannpolynomial ˆ 2 (2 2 ) ˆ +(1 ).Wecancheckthat therootsofthispolynomialwillbecomplexif0 < p 2+2 p 1 ,andthatinthiscase therootssatisfy j ˆ j 2 = 1 4 (2 2 ) 2 +4(1 ) (2 2 ) (2.7) =1 < 1(2.8) whichshowsboththestabilityanddissipativenatureofthesecondorderscheme. 34 FourthOrderScheme WeobtaintheVonNeumannpolynomial ˆ 2 (2 2 2 + 4 = 12) ˆ +(1 ).Wecancheck thattherootsofthispolynomialwillbecomplexif0 < s 12 2 q 1 1 6 1+ p 1+ , andthatinthiscasetherootssatisfy j ˆ j 2 = 1 4 (2 2 2 + 4 = 12) 2 +4(1 ) (2 2 2 + 4 = 12) (2.9) =1 < 1(2.10) whichshowsboththestabilityanddissipativenatureofthefourthorderscheme. Wenotethatineachcase,themaximumallowedvalueof isslightlysmallerthanwhat isallowedbythecorrespondingschemewithoutdissipation.Amoredetailedanalysisshows thattheedampingrateis k 2 k 2 + 2 2 forthesecondorderscheme,and k 2 k 2 + 2 3 forthefourthorderscheme,meaningthathighfrequenciesaremorerapidlydampedthan lowfrequencies. 2.6.2Dissipationin2D Againmotheschemesgivenin[11],wehavethesecondorderschemewithdissipation, u n +1 2 u n + u n 1 = 2 C [ u n ]+ C 2 [ u n 1 ] ; (2.11) andthefourthorderschemewithdissipation, 35 u n +1 2 u n + u n 1 = 2 C [ u n ] 2 DC [ u n ] 4 12 C 2 [ u n ]+ C 3 [ u n 1 ] : (2.12) wherenow D =1 L 1 x L 1 y and C = L 1 y D x + L 1 x D y .Forfurtherdetailsonthese operators,see[11].Numericalexperimentsindicatethatthe2Dschemeswithar dissipationareindeedunconditionallystablewiththesamemaximumvalueof aswiththe 1Dschemes. 2.7DerivationofaFourth-OrderDissipativeNewmark Scheme AsanalternativetotheartidissipationtermsforcenteredschemesgiveninSection 2.6,wenowgivethederivationofafourth-orderaccurate,dissipative,and(apparently) unconditionallystablemethodthatwecallaNewmarkscheme,duetotheintroductionof anauxiliaryvariable v = u t . Considerasmoothsolution u ( x;t )to(2.1).WecanformtheTaylorexpansion u ( x;t + t )= p 1 X k =0 @ k u @t k ( x;t ) t k k ! + @ p u @t p ( x; ( x;t )) t p p ! Thentoforma p -thorderaccurateapproximationto u ( x;t + t ),wemustform( p k )-th orderapproximationsto @ k u @t k ( x;t ),for k =0 ;:::;p 1. Alternatively,wecouldkeeptheoddtimederivativesin(2.7),andconsidertheTaylor 36 expansionof v = u t ,usingthefactthat v solvethewaveequation v tt c 2 v = f t ( x;t ): v ( x;t + t )= 1 X k =0 @ k v @t k ( x;t ) t k k ! = 1 X k =0 @ 2 k v @t 2 k ( x;t ) t 2 k (2 k )! + + 1 X k =1 @ 2 k u @t 2 k ( x;t ) t 2 k 1 (2 k 1)! = v ( x;t )+ 1 X k =1 2 4 c 2 k k v ( x;t )+ k X j =1 c 2( k j ) k j @ 2 j 1 @t 2 j 1 f ( x;t ) 3 5 t 2 k (2 k )! + + 1 X k =1 2 4 c 2 k k u ( x;t )+ k X j =1 c 2( k j ) k j @ 2( j 1) @t 2( j 1) f ( x;t ) 3 5 t 2 k 1 (2 k 1)! : WritingoutthefewtermsoftheTaylorexpansionsexplicitly: u ( x;t + t )= u ( x;t )+ tv ( x;t )+ t 2 2 c 2 u + f ( x;t )+ + t 3 6 ( c 2 v + f t )( x;t )+ O t 4 ) v ( x;t + t )= v ( x;t )+ t c 2 u + f ( x;t )+ t 2 2 ( c 2 v + f t )( x;t )+ + t 3 6 c 4 2 u + c 2 f + c 2 f tt ( x;t )+ O t 4 ) Inourtargetapplication(particlesimulationofplasmas),thesourcetermrepresentsthe contributionofparticles,andmustbehandledwithsomecare.Fornow,weignorethesource 37 term: u ( x;t + t )= u ( x;t )+ tv ( x;t )+ t 2 2 c 2 u ( x;t )+ t 3 6 c 2 v ( x;t )+ O t 4 ) v ( x;t + t )= v ( x;t )+ tc 2 u ( x;t )+ t 2 2 c 2 v ( x;t )+ + t 3 6 c 4 2 u ( x;t )+ O t 4 ) InordertoobtainanA-stablescheme,wereplacetheLaplacianoperatoranditspow- ers,whichareunbounded,withboundedapproximatingoperators.Weconsidertheone- dimensionalcase,withhigher-dimensionstobehandledthroughdimensionalsplitting.We usetheidentityfrom[11]: @ xx 2 m =( 1) m 1 X p = m p 1 m 1 D p where D [ u ]( x )= u ( x ) 2 Z 1 e j x y j u ( y ) dy Wewriteoutthefewtermsoftheapproximationsweneed: 38 @ xx 2 = D 2 + @ xx 2 2 = D 2 +2 D 3 + Inserting: u ( x;t + t )= u ( x;t )+ tv ( x;t )+ t 2 2 c 2 2 ( D 2 ) u ( x;t )+ + t 3 6 c 2 2 ( D 2 ) v ( x;t )+ O t 4 ) v ( x;t + t )= v ( x;t )+ tc 2 2 ( D 2 ) u ( x;t )+ t 2 2 c 2 2 ( D 2 ) v ( x;t )+ + t 3 6 c 4 4 ( D 2 +2 D 3 ) u ( x;t )+ O t 4 ) Choosing = =c t : u ( x;t + t )= u ( x;t )+ tv ( x;t )+ 2 2 ( D 2 ) u ( x;t )+ + 2 6 ( D 2 ) v ( x;t )+ O t 4 ) v ( x;t + t )= v ( x;t )+ 2 t ( D 2 ) u ( x;t )+ 2 2 ( D 2 ) v ( x;t )+ + 4 t ( D 2 +2 D 3 ) u ( x;t )+ O t 4 ) Numericalexperimentsindicatethatthefollowingscheme,includingselecthigher-order terms,maintainsunconditionalstabilityfor0 < 1,withasmallamountofdissipation: 39 u n +1 = u n + tv n 2 2 D u n + D 2 u n (2.13) 2 t 6 D v n + D 3 v n + 4 24 D 2 u n v n +1 = v n 2 t D u n + D 2 u n 2 2 D v n + D 2 v n + D 3 v n + + 4 t D 2 u n +2 D 3 u n + 4 24 D 2 v n Thisschemepossessescertaindisadvantagescomparedtothecenteredschemeswith dissipation.Namely,itrequiresgreaterstorageandcomputationalrequirements (requiringalargernumberofconvolutions),andthelevelofdissipationisnotobviously tunable.Ontheotherhand,itmaybepossibletofurthermodifythecotsofthe higher-ordertermstoallowfortuningofthelevelofdissipation,andfurther,inelectromag- neticPICsimulations,weneedtocomputethetimederivative A t inordertocalculatethe electricwhichwouldbecomputedviainthecaseofcenteredschemes, butwouldbegivenautomaticallybytheNewmarkmethodtofourth-orderaccuracy. 2.8EmbbededBoundaryMethodforNeumannBound- aryConditions InimplementingNeumannboundaryconditionsforboundarygeometriesconformingtogrid lines,suchasarectangulardomain,wecandirectlyimposeatwo-pointboundarycorrection. Onewaytoextendthismethodtoageneralpolygonaldomainwouldbetousemultiple oversetgrids,eachalignedwithaboundarysegment,whichcommunicatewiththeinterior gridthroughinterpolationonaghostcellregion,thoughwedonotpursuethatapproachin 40 thiswork. Forcurvedboundaries,analternativeapproachisthatofanembeddedboundarymethod, whichinvolvesdeterminingtheDirichletvaluesattheendpointsofeach x -and y -sweeplines thatresultintheapproximatesatisfactionoftheNeumannboundarycondition(in constructinganapproximateNeumann-to-Dirichletmap).Wepresenttheimplementation ofanembeddedboundarymethodforNeumannboundaryconditionsfortheimplicitwave solveronacurvedboundarygeometry.Theapproachtakenherefollowstheworkin[28], whichproposesanembeddedboundarymethodforNeumannboundaryconditionswitha methodforthewaveequation.Theanalysisinthatworksuggeststhat,on thecontinuouslevel,themoequationsandboundaryconditionsresultingfromtypical truncationerrortermspossessunstableboundarylayersolutions,sothattheadditionof adissipativetermisnecessarytoachieveastablemethod.Thisisconsistentwithour experienceintheimplementationdescribedhere,withtheembeddedboundarymethod becomingunstablewhenappliedtothenon-dissipativedispersivesolver,butremainingstable fortheesolver,whichisdissipative. Inthefollowing,webriedescribethetwo-pointboundarycorrectionmethodfora grid-alignedrectangularboundary.Wethendescribetheembeddedboundarymethodfor a1Dproblem.Thismethodrequiresaniterativeprocedure,whichweshow,inthesetting ofthe1Dproblem,tobeaconvergentcontractionmappingwitharateofconvergencethat dependsontheCFLnumber.Wedescribetheimplementationoftheembeddedboundary methodin2D,andgivenumericalresults. 41 2.8.1DescriptionoftheTwo-PointBoundaryCorrectionMethod Inarectangulardomainwheretheboundariesconformtogridlines,itisstraightforwardto imposethetwo-pointboundarycorrectiontermsinaline-by-linefashion,sinceinthiscase, thegridlinesarenotcoupledthroughthenormalderivative.Asthisisasimpleextension ofthe1Dboundarycorrectionalgorithm,wedonotelaboratefurther. 2.8.2DescriptionoftheEmbeddedBoundaryMethodandProof ofConvergenceoftheIterativeSolutionin1D Weconsiderthesituationofaone-dimensionaldomain f x B 0and II = ˘ 2 G ˘ 2 I ˘ 2 II ˘ 2 I < 0are O (1),weonlyneed supplysecond-orderaccurateapproximationsto u I and u II tomaintainoverallsecond-order accuracy.(Thecotswouldbe O (1 = x )inthecaseofnonhomogeneousNeumann boundaryconditions,whichwouldrequirethird-orderaccurateapproximationsto u I and u II .Forsimplicity,weconsideronlythecaseofhomogeneousNeumannboundaryconditions inthepresentwork.)Suchapproximationsmaybeobtainedthroughlinearinterpolation, giving 43 u I = ˙ I u m +(1 ˙ I ) u m +1 u II = ˙ II u n +(1 ˙ II ) u n +1 where ˙ I = x m +1 x I x , ˙ II = x n +1 x II x ,and u j = u ( x j )arethevaluesofthefunction attheuniformgridpointsfor j = m;m +1 ;n;n +1. x 0 x G x B x 1 x m x m +1 x I x n x n +1 x II Figure2.1:Boundarygeometryin1D. Hence,todeterminetheghostpointvalue u G thatleadstoasolutionconsistentwith homogeneousNeumannboundaryconditions,weshouldsolvethefollowingsystemofequa- tions. u G = I ( ˙ I u m +(1 ˙ I ) u m +1 )+ II ( ˙ II u n +(1 ˙ II ) u n +1 ) u m = I m +( u G I G ) e ( x m x 0 ) u m +1 = I m +1 +( u G I G ) e ( x m +1 x 0 ) u n = I n +( u G I G ) e ( x n x 0 ) u n +1 = I n +1 +( u G I G ) e ( x n +1 x 0 ) with I , II , ˙ I and ˙ II asabove,andwhere I j = I ( x j )istheconvolution integralevaluatedatuniformgridpointsfor j = m;m +1 ;n;n +1.Thissystemcanbe 44 solvedformallyfor u G ,giving u G = h I ˙ I ( I m I G e ( x m x G ) )+(1 ˙ I )( I m +1 I G e ( x m +1 x G ) ) + II ˙ II ( I n I G e ( x n x G ) )+(1 ˙ II )( I n +1 I G e ( x n +1 x G ) ) h 1 I ˙ I e ( x m x G ) +(1 ˙ I ) e ( x m +1 x G ) II ˙ II e ( x n x G ) +(1 ˙ II ) e ( x n +1 x G ) Toshowthatthissolutionformulaiswwearguethat 0 0, II < 0,0 ˙ I 1,and 0 ˙ II 1,weobtain K = I d m [ ˙ I +(1 ˙ I ) d ]+ II d n [ ˙ II +(1 ˙ II ) d ] I d m + II d n +1 = d m ˘ 2 II d n +1 ˘ 2 I + ˘ 2 G ( d n +1 d m ) ˘ 2 II ˘ 2 I d m ˘ 2 II d n +1 ˘ 2 I ˘ 2 II ˘ 2 I = s 2 I d m d n +1 s 2 I s 2 I s 2 I = 4 d m d n +1 3 45 Now,since x< s I < (3 = x ,wecanseethatitisthecasethateither m =1and n =2or3,orthat m =2and n =3.Itisthenamatterofsomesimplecalculustocheck thatthatthefunctions f m;n ( x )=(4 x m x n +1 ) = 3satisfy f m;n ( x ) < 1for0 0 Inthetwo-dimensionalcase,theline-by-linesolutionmethodcouplestheghostpoint values,andageneralexplicitsolutionformulaisimpossibletowritedown.Inprinciple,one maywriteoutanddirectlysolvealinearsystemtoobtaintheghostpointvalues.Instead,we proposeaniterativesolutionmethodthatavoidstheformationofamatrix.Wenowdescribe thisiterativesolutionmethodandproveitsconvergenceinthecontextoftheone-dimensional problemdescribedabove. Supposewehavetheconvolutionintegralevaluatedatthegridpoints, I j ,anda k -th 46 iteratefortheghostpointvalue, u k G .Thenwemayobtainthenextiteratebytheformulas u k +1 m = I m +( u k G I G ) e ( x m x 0 ) u k +1 m +1 = I m +1 +( u k G I G ) e ( x m +1 x 0 ) u k +1 n = I n +( u k G I G ) e ( x n x 0 ) u k +1 n +1 = I n +1 +( u k G I G ) e ( x n +1 x 0 ) u k +1 G = I ( ˙ I u k +1 m +(1 ˙ I ) u k +1 m +1 )+ II ( ˙ II u k +1 n +(1 ˙ II ) u k +1 n +1 ) wherequantitiesareasabove.Now,toprovetheconvergenceoftheinteration,we showitiscontractive.Takingtheoftwoiterates,wehave j u k +1 G u k G j = j I ( ˙ I ( u k +1 m u k m )+(1 ˙ I )( u k +1 m +1 u k m +1 ))+ II ( ˙ II ( u k +1 n u k n )+(1 ˙ II )( u k +1 n +1 u k n +1 )) j = j I ˙ I ( u k G u k 1 G ) e ( x m x G ) +(1 ˙ I )( u k G u k 1 G ) e ( x m +1 x G ) + + II ˙ II ( u k G u k 1 G ) e ( x n x G ) +(1 ˙ II )( u k G u k 1 G ) e ( x n +1 x G ) j K j u k G u k 1 G j where0 0and II = ˘ 2 G ˘ 2 I ˘ 2 II ˘ 2 I < 0.Inthe2D case,weapproximationsto u I and u II throughbilinearinterpolation.Thisisincontrast 48 to[28],whotheintersectionofthenormalwithgridlines,theninterpolatealongthe gridlines.Wehavealsoimplementedasecond-orderaccurateversionofthisapproachand comparedtothebilinearinterpolationschemeproposedhere.Wehavefoundthatthetwo schemesbehavesimilarly,howeverthebilinearinterpolationschemeisslightlymoreaccurate andsimplertocode,notrequiringtohandleseparatecasesofintersectionwithhorizontal andverticalgridlines.Thebilinearinterpolationschemeisstandard,butwegiveithere forcompleteness.Iftheinterpolationpoint u I liesinacellwithcorners( x i ;y j ),( x i +1 ;y j ), ( x i +1 ;y j +1 )and( x i ;y j +1 ),thenwehavethefollowingapproximationfor u I : u I = w 1 u i;j + w 2 u i +1 ;j + w 3 u i +1 ;j +1 + w 4 u i;j +1 where w 1 = ( x i +1 x I )( y j +1 y I ) x y , w 2 = ( x I x i )( y j +1 y I ) x y , w 3 = ( x I x i )( y I y j ) x y and w 4 = ( x i +1 x I )( y I y j ) x y .Withthisinterpolationschemeestablished,wenowoutlinethe algorithmforthe2Ddimensionally-splitwavesolver. Theaboveinterpolationprocedureappliesregardlessofthevarietyofthewavesolver thatitisusedwith,providedthatthewavesolverhastdissipationtomaintain stability.Wenowdescribetherestoftheembeddedboundaryalgorithminthecontextof theewavesolver,thoughititmaybesimilarlyappliedtothethedispersiveschemes withdissipationorthedissipativeNewmarkschemedescribedabove.Inanalogy totheiterationpresentedinthe1Dcase,the2Diterativealgorithmproceedsbyusingthe interpolationschemetoprovidevaluesattheghostpoints,whichinturnprovidenewvalues fortheboundarycorrectioncots,whicharethenusedtoupdatethevaluesatthe interiorgridpoints,comprisingonefulliteration.Itshouldbenotedthatnotallinterior 49 gridpointsneedbeupdatedintheiteration,onlythoseneartheboundarythatliewithin theboundaryinterpolationstencils. Usingvaluesfromprevioustimesteps,theinitialguessfortheinteriorgridpointsin theboundaryinterpolationstencilsmaybegivenbyextrapolationintime,as u n +1 ; 0 = 2 u n u n 1 (linearextrapolation)or u n +1 ; 0 =3 u n 3 u n 1 + u n 2 (quadraticextrapolation). Eitherextrapolatedinitialguessprovidesamodestreductioninthenumberofiterations requiredversusazeroinitialguess,withonlyaslightfurtherreductioninthenumberof iterationsgoingfromlineartoquadraticextrapolation.Aneestoppingcriterionfor iterationis j u n +1 ;l +1 u n +1 ;l j 1 < ,where issomechosentolerance,whichmaybechosen tobequitesmall,astheiterationisapointinteration.Inthenumericalexample,we chooseatoleranceof10 15 ,andweachieveconvergenceinlessthan40iterationsataCFL numberof2. Inapplyingtheeversionofthewavesolver,weassumewehaveprevioustimesteps u n , u n 1 and u n 2 .WehavetosolvethemoHelmholtzequationwithhomogeneous Neumannboundaryconditions, 1 1 2 r 2 u n +1 = 1 2 5 u n 4 u n 1 + u n 2 in @u @n =0on= @ where = p 2 c t .Weapplydimensionalsplittingto 1 1 2 r 2 u n +1 = 1 1 2 @ xx 1 1 2 @ yy u n +1 + O 1 4 50 sowed w = 1 1 2 @ yy u ,andnotingthat w = u + O ( c t ) 2 sothat @w @n = @u @n + O ( c t ) 2 ,weobtainthefollowingapproximatesystem 1 1 2 @ xx w n +1 = 1 2 5 u n 4 u n 1 + u n 2 in @w @n n +1 =0on= @ 1 1 2 @ yy u n +1 = w n +1 in @u @n n +1 =0on= @ WenowsupposeourdomainisembeddedinauniformCartesiangrid,withhorizontal gridlinescorrespondingto y = y k ,1 k N y andverticalgridlinescorrespondingto x = x j ,1 j N x .Theembeddedboundaryalgorithmwillbeappliedwhencalculating theintermediatevariable w n +1 inhorizontallinesweepsaswellasthesolutionvariable u n +1 inverticallinesweeps.Theiterationsforthesetwovariablesareseparate;theiterative procedureisappliedto w toconvergence,andthenthisvalueof w isusedtocompute u ,and theiterativeprocedureisappliedto u toconvergence.However,ineachiteration,thegrid linesarecoupledthroughtheinterpolationscheme,sothatallgridlinesmustbeiterated together.Theoveralliterativealgorithmisdescribedin1,withdetailsspforthe iterationon w .Theiterationon u isverysimilar,andsoweomitthedetails. 51 Algorithm1 ApplicationofNeumannBoundaryConditionsin2dimensions 1. (Initializationofghostpoints) Performtheinterpolationschemedescribedabove toobtainthevaluesof u n , u n 1 ,and u n 2 attheghostpoints,whicharetheendpoints ofthehorizontalandverticalgridlines. 2. (Evaluationofparticularsolution) Foreachhorizontalline y = y k ,for1 k N y , withghost(end)points x = a k and b k theparticularsolution w n +1 p;k forthe intermediatevariable w n +1 k ( x )byevaluatingthediscreteconvolutionoperator w n +1 p;k ( x j )= 4 Z b k a k [5 u n 4 u n 1 + u n 2 ]( x 0 ;y k ) e j x j x 0 j dx 0 foreachgridpoint x j inthehorizontalline,includingtheghostpoints. 3. (Boundarycorrectioninitialization) Foreachhorizontalline y = y k ,settheinitial guessfortheintermediatevariableviaextrapolation, w n +1 ; 0 k =3 w n k 3 w n 1 k + w n 2 k , ontheinteriorpointswithintheboundaryinterpolationstencil. 4. (Boundarycorrectioniteration) Foreachhorizontalline y = y k ,performthe interpolationschemeusingtheinteriorvaluesof w n +1 ;l p;k totheghostpointvalues. Usingtheseghostpointvalues,applytheboundarycorrectiononeachlinetoobtain theupdatedintermediatevariable, w n +1 ;l +1 k ( x j )= w n +1 p;k ( x j )+ A k e ( x j a k ) + B k e ( b k x j ) forthevaluesof x j lyingwithintheboundaryinterpolationsten- cil,where A k = w n +1 ;l k ( a k ) w n +1 p;k ( a k ) k w n +1 ;l k ( b k ) w n +1 p;k ( b k ) 1 2 k , B k = w n +1 ;l k ( b k ) w n +1 p;k ( b k ) k w n +1 ;l k ( a k ) w n +1 p;k ( a k ) 1 2 k , k = e ( b k a k ) .Checkfor convergence,andifconverged,storetheintermediatevariable w n +1 . 5. Repeatthisprocessfortheverticallinesweeps,usingtheintermediatevariable w n +1 to calculatetheparticularsolutionfor u n +1 ,thenapplytheboundaycorrectioninteration. 52 MoHelmholtzEquation DispersiveScheme, = 2 c t : 1 2 r 2 +1 u n +1 +2 u n + u n 1 =4 u n + 4 2 f ( x;t n ) eScheme, = p 2 c t : 1 2 r 2 +1 u n +1 = 1 2 5 u n 4 u n 1 + u n 2 + 1 2 f ( x;t n +1 ) DimensionallySplit 1 2 r 2 +1 u = f ) MoHelmholtzEquation 1 2 @ 2 @x 2 +1 1 2 @ 2 @y 2 +1 u = f ) 1 2 @ 2 @x 2 +1 w = f , 1 2 @ 2 @y 2 +1 u = w 1DIntegralSolution 1 2 d 2 dx 2 +1 u = f on( a;b ) ) u ( x )= 2 R b a f ( x 0 ) e j x x 0 j dx 0 + Ae ( x a ) + Be ( b x ) = I [ f ]( x )+ Ae ( x a ) + Be ( b x ) 1DBCCorrectionCots Dirichlet: A = ( u a I a ) ( u b I b ) 1 2 , B = ( u b I b ) ( u a I a ) 1 2 u ( a )= u a , u ( b )= u b Neumann: A = ( v b + b ) ( v a a ) (1 2 ) , B = ( v b + b ) ( v a a ) (1 2 ) u 0 ( a )= v a , u 0 ( b )= v b Periodic: A = I b 1 , B = I a 1 u ( a )= u ( b ), u 0 ( a )= u 0 ( b ) I a = I [ f ]( a ), I b = I [ f ]( b ), = e ( b a ) FastConvolutionAlgorithm a = x 0 > > < > > > : cos (2 m +1) ˇx L x cos (2 n +1) ˇy L y Dirichletcase sin (2 m +1) ˇx L x sin (2 n +1) ˇy L y Neumanncase and u t ( x;y; 0)=0 for( x;y ) 2 m and n integers.Exactsolutionsarewell-knownineachcase.Theresults ofementstudiesarelistedinTables2.2and2.3.Theerroristhemaximumdiscrete L 2 error(computedagainsttheexactsolution)overalltimesteps. 54 ( x G ;y G ) ( x I ;y I ) ( x II ;y II ) ( x B ;y B ) Figure2.2:Boundarygeometryin2D. x y y =0 x =0 Figure2.3:Rectangularcavitywithdomaindecomposition. CFL0.5 CFL2 CFL10 x L 2 error L 2 order L 2 error L 2 order L 2 error L 2 order 1 = 40 8 : 4801 10 4 5 : 8919 10 3 1 : 0692 10 1 1 = 80 1 : 9018 10 4 2 : 15671 1 : 4578 10 3 2 : 01499 3 : 3924 10 2 1 : 65620 1 = 160 4 : 4811 10 5 2 : 08544 3 : 6061 10 4 2 : 01523 8 : 8683 10 3 1 : 93555 1 = 320 1 : 0861 10 5 2 : 04468 8 : 9550 10 5 2 : 00965 2 : 2080 10 3 2 : 00592 1 = 640 2 : 6726 10 6 2 : 02284 2 : 2300 10 5 2 : 00535 5 : 4638 10 4 2 : 01476 Table2.2:tstudyforrectangularcavitywithDirichletBCusingdomaindecom- position.Here, c =1, m = n =0,and L x = L y =1. 55 CFL0.5 CFL2 CFL10 x L 2 error L 2 order L 2 error L 2 order L 2 error L 2 order 1 = 40 8 : 8928 10 4 6 : 1862 10 3 1 : 1147 10 1 1 = 80 1 : 9486 10 4 2 : 1902 1 : 4942 10 3 2 : 0497 3 : 4771 10 2 1 : 6807 1 = 160 4 : 5367 10 5 2 : 1027 3 : 6511 10 4 2 : 0329 8 : 9791 10 3 1 : 9532 1 = 320 1 : 0929 10 5 2 : 0535 9 : 0110 10 5 2 : 0185 2 : 2217 10 3 2 : 0148 1 = 640 2 : 6809 10 6 2 : 0273 2 : 2370 10 5 2 : 0098 5 : 4800 10 4 2 : 0192 Table2.3:tstudyforrectangularcavitywithNeumannBCusingdomaindecom- position.Here, c =1, m = n =0,and L x = L y =1. 2.9.2DoubleCircleCavity Inthisexample,wesolvethewaveequationwithhomogeneousDirichletboundaryconditions ona2Ddomainwhichis,asinFigure2.4,theunionoftwooverlappingdisks,withcenters P 1 =( ; 0)and P 2 =( ; 0),respectively,andeachwithradius R : = f ( x;y ): j ( x;y ) P 1 j > > > > > > > < > > > > > > > > : cos 6 ˇ 2 j ( x;y ) P 1 j 0 : 8 2 j ( x;y ) P 1 j < 0 : 8 cos 6 ˇ 2 j ( x;y ) P 2 j 0 : 8 2 j ( x;y ) P 2 j < 0 : 8 0otherwise and u t ( x;y; 0)=0 for( x;y ) 2 SelectedsnapshotsoftheevolutionaregiveninFigure2.5,andtheresults ofatstudyaregiveninTable2.4.Thediscrete L 2 errorwascomputedagainst awnumericalreferencesolution x =2 : 1875 10 4 );theerrordisplayedinthe tableisthemaximumovertimestepswith t 2 [0 : 28 ; 0 : 29].Forthisexample, R =0 : 3, =0 : 2, c =1,andtheCFLis2. x y t L 2 error L 2 order 7 : 0000 10 3 4 : 3333 10 3 8 : 6667 10 3 6 : 1437 10 3 3 : 5000 10 3 2 : 1667 10 3 4 : 3333 10 3 1 : 6829 10 3 1 : 8681 1 : 7500 10 3 1 : 0833 10 3 2 : 1667 10 3 4 : 3595 10 4 1 : 9488 8 : 7500 10 4 5 : 4167 10 4 1 : 0833 10 3 1 : 0515 10 4 2 : 0517 Table2.4:tstudyforthedoublecirclecavitywithDirichletBC.Forthenumerical referencesolution, x =2 : 1875 10 4 , y =1 : 3542 10 4 ,and t =2 : 7083 10 4 . 2.9.3PeriodicSlitGrating Inthisexample,weapplyourmethodtomodelanperiodicgrating underanincidentplanewave.gratingsareperiodicstructuresusedinoptics toseparatetwavelengthsoflight,muchlikeaprism.Thehighresolutionthatcan 58 beachievedwithgratingsmakesthemusefulinspectroscopy,forexample,in thedeterminationofatomicandmolecularspectra.Ournumericalexperiment,depicted inFigure2.6,demonstratestheuseofourmethodwithmultipleboundaryconditionsand nontrivialgeometryinasinglesimulationtocapturecomplexwavephenomena. Anidealizedslitractiongratingconsistsofascreenofvanishingthickness, withopenslitsofaperturewidth a ,spaceddistance d apart,measuredfromtheendofone slittothebeginningofthenext(thatis,theperiodicityofthegratingis d ).Weimposean incidentplanewaveoftheform u inc ( x;y;t )=cos( !t + ky ),where k =2 ˇ=a and ! = k=c , where c isthewavespeed.PeriodicBCsat x = d= 2(determiningtheperiodicityofthe grating),andhomogeneousDirichletBCsareimposedatthescreen.Wealsotestw boundaryconditionsinmultipledimensions,whichareimposedat y = L y = 2. x y a d u inc wBC wBC PeriodicBC PeriodicBC Figure2.6:Periodicslitgratinggeometry InFigure2.7,weobservethetimeevolutionoftheincidentplanewavepassingthrough theaperture,andtheresultinginterferencepatternsasthewavepropagatesacross theperiodicboundaries.Thewboundaryconditionsallowthewavestopropagate outsidethedomain.WhilearigorousanalysisoftheofourwBCsisthe subjectoffuturework,theresultslookquitereasonable,asnospuriousareseen attheboundaries. 59 2.9.4BesselModewithNeumannBoundaryConditions Herewepresentanumericalexampleoftheembeddedboundarymethodforhomogeneous NeumannboundaryconditionsgiveninSection2.8.Weapplythemethodtoacircular domain,forwhichanalyticalsolutionsexist.Weconsideraradially-symmetricBesselmode withhomogeneousNeumannboundaryconditions,withananalyticsolutiongivenby u ( r;t )= J 0 Z 0 r R cos Z 0 ct R ; (2.14) where J 0 istheBesselfunctionofthekindoforder0, r = p x 2 + y 2 , R isthe radiusofthedomain,and Z 0 ˇ 3 : 8317isthesmallestnonzerorootof J 0 0 (sothat @u @ n ( R;t )= @u @r ( R;t )= Z 0 R J 0 0 ( Z 0 )cos Z 0 ct R =0).Inthisexample,wetakeradius R = ˇ= 2andwave speed c =1.AnexampleoftheembeddedboundarygridusedisgiveninFigure2.8.We performatstudywithaCFLnumberof2,withtheresultsinFigure2.9 indicatingtheexpectedsecond-orderconvergence.Wesettheiterationtoleranceto10 15 , andweseeconvergenceoftheboundarycorrectioniterationinfewerthan40iterations.We notesomeoscillationofthe L 1 erroraboutthelinegivingsecond-orderaccuracy,which webelievetobeduetothegridpointsmovingwithrespecttotheboundarythroughthe t,causingsomevariationinthemaximumerror. 60 (a) t =0 : 31 (b) t =0 : 51 (c) t =1 : 01 (d) t =2 : 01 Figure2.7:Evolutionoftheslitgratingproblem,withaperturewidth a =0 : 1, gratingperiodicity L y = d =1,andwavespeed c =1.TheCFLisat2. Figure2.8:Anexampleoftheembeddedboundarygridused.Theredcircledexteriorgrid pointsaretheendpointswhereavalueistobecalculatedviatheinterpolationprocedure. Theredcrossesarethepointswherevaluesareimposedonquadraticboundaryinterpolant alongthenormaldirection(reddashedline).Valuesforthebilinearinterpolantsaresupplied fromthegreencircledinteriorgridpoints. 61 Figure2.9:tstudyfortheBesselmodeinacirculardomainwithCFL numberof2.Usingquadraticboundaryinterpolantwithbilinearinteriorinterpolant. 62 Chapter3 Descriptionoftheparticle-in-cell method 3.1ParticleWeightingScheme InourPICmethods,thechargeandcurrentdensitiesarerepresentedasthesumofparticle shapefunctions: ˆ ( x ;t )= N p X i =1 q i S x x p;i ( t ) x (3.1) J ( x ;t )= N p X i =1 q i v p;i S x x p;i ( t ) x (3.2) where N p isthenumberofparticles, x p;i ( t )and v p;i ( t )arethepositionandvelocity, respectively,ofparticle i attime t ,and S ( x )isaparticleshapefunction.Itshouldbe emphasizedthatthesearenotphysicalparticles,butrathermacro-orsuperparticlesthat representadiscretizationofthePDF[4].Weanalyticallyevaluatetheparticleconvolution integralscorrespondingtothesesourcetermswiththealgorithmsdescribedbelow. 63 3.2MethodforControllingDivergenceError ItiswellknownthatelectromagneticPICmethodswhichdonotsatisfyadiscreteformof Gauss'lawthroughtheirsolversandchargeandcurrentweightingschemeswill severenumericalerrorsrelatedtochargeconservation[37].Weseektodevelopstaggeredgrid approachestocomputingpotentialsandldsthatsatisfydiscreteanaloguesofGauss'law r E = 0 andtheidentity r B =0.Thedivergence-freeconditionfor B willbeeasily ingeneral,asthemagneticwillbecalculatedasacurlofthe vectorpotential.InordertonumericallyenforceGauss'law,weseektoperformanelliptic divergencecorrection.Futureworkwillconsiderthealternativeofhyperbolicdivergence cleaning[41]. Wenowgivethemathematicalunderpinningofourellipticdivergencecorrectiontech- nique.Theelectriciscalculatedas E = @ A @t .ThenGauss'lawmayberewritten as: 0 = r E (3.3) = r ( @ A @t )(3.4) = @ ( r A ) @t : (3.5) Thus,thescalarpotentialthePoissonequation = 0 + @ ( r A ) @t (3.6) OurmethodisbasedontheobservationthatifthisPoissonequationissuitablydis- 64 cretizedandsolvedonastaggeredgridtoprovidethescalarpotentialusedincalculating theelectricthentheelectricwillautomaticallysatisfyadiscreteformofGauss'law. Whiletheexactformofthestaggeredgridwilldependonwhichcomponentsofthecurrent densityandtheelectricandmagneticeldsareretainedinagivenmodel,ourmethodguar- anteesexactdiscretedivergencerelationsindependentlyofthechargeandcurrentweighting schemesused,ofthenatureofthesolverusedfor A ,andofthegaugeconditionsp. Detailsaboutthespstaggeredgridsusedaregivenin3.6.2,alongwithproofsofthe exactdiscretedivergencerelations.Itshouldbenotedthatthisprocedureisananalogueof theellipticdivergencecorrectiontechniquespresentedin[32,37],andalsobearssimilarityto themethodgivenin[40]toenforcethedivergencerelationswhenevolvingthepotentialsin theCoulombgauge.Toimposeoutwboundaryconditionsonthescalarpotential,wesolve anauxiliarywaveequationforusingourwavesolverwithwboundaryconditions, whichthensuppliestheboundaryvaluesforinthePoissonsolve. 3.3ParticleEquationsofMotion InourPICmethods,theapproximationoftheevolutionoftheVlasovequationamountsto theintegrationoftheequationsofmotionoftheparticles: d x p;i dt = v p;i ( t )(3.7) d v p;i dt = a p;i ( t )(3.8) where a p;i ( t )istheaccelerationofparticle i attime t .Toevolvetheparticleequationsof 65 motion,weobtainparticleaccelerationsthroughtheusualinterpolationoffromgrid pointstoparticlelocations[4],andweusestandardnumericalmethods,suchastheexplicit leapfrogmethodforelectrostaticproblemsandtheforwardEulermethodforelectromagnetic problems(theimplicitsourcetermintheewavesolvercausessomeyin achievinghigher-orderaccuracyintheintegrationoftheparticleequationsofmotion-for thisreasonthedispersiveschemeandhigher-ordercenteredschemesgivenin[11]willbe pursuedinfuturework-buttheaccuracyachievedistforthenumericalresults inthiswork).Fieldsarecalculatedasofpotentialsonthegrid.Asthese aspectsarestandardandnotthefocusofthepresentwork,wedonotelaboratefurther. 3.4FastConvolutionAlgorithmForParticleSources 3.4.1FastConvolutionAlgorithmin1D Wenowdescribethealgorithmusedforthefastexactevaluationoftheconvolutionof chargeandcurrentdensitysourcetermsduetoparticles.Ithastwomainsteps.Thereis alocaldepositstepandthenarecursivesweepstep.Thisbasicstructureisthesameinall dimensions.Forwedescribetheapplicationofthealgorithmtolinearparticle shapesinoneandtwodimensions.However,itmaybegeneralizedtoanyseparableparticle shapeswithcompactsupportinanydimension.Notethatthisincludesmanywidelyused particleshapesinPICalgorithms,namelytypicalspline-basedparticleshapesand(suitably )Gaussianparticleshapes.Forthecaseofthechargedensityintegral,theparticle shapefunction S p belowisreplacedby q p S p anditscontributionsummedtothecharge densityintegral,andforthecaseofthecurrentdensityintegral, S p isreplacedby v p q p S p anditscontributionsummedtothecurrentdensityintegral. 66 3.4.1.1LocalDepositStepin1D Considerparticle p locatedinthecell[ x m ;x m +1 ]inauniformgridwithcelllength x . Let S p ( x )betheshapefunctionoftheparticle.Assumethatthesupportof S p haslength 2 r x forsomeinteger r .Thelocaldepositstepthenconsistsofanalyticallyevaluatingthe integrals J L;p j +1 = Z x j +1 x j S p ( x 0 ) e ( x j +1 x 0 ) dx 0 (3.9) J R;p j = Z x j +1 x j S p ( x 0 ) e ( x 0 x j ) dx 0 (3.10) for j = m r;:::;m;:::;m + r andforeachparticle p ,andsummingtheirvaluesonthe grid. Forlinearparticleshapes(correspondingto r =1),wehave S p ( x )= S ( x x p x ),where S ( x )= 8 > > > < > > > : 1 j x jj x j < 1 ; 0 j x j 1 (3.11) Let a =( x p x m ) = x ,where x p isthelocationoftheparticle.Forsimplicity,let x p =0 and x =1.Forlinearparticleshapes,wethenhavethesituationdisplayedinFigure3.1. Thedesiredintegralsaretheneasilyevaluatedforlinearparticleshapes: 67 J L;p m = Z a 1 (1+ x 0 ) e ( a x 0 ) dx 0 (3.12) =(((1 a ) 1)+ e a e ) (3.13) J L;p m +1 = Z 1 a a (1 j x 0 j ) e (1 a x 0 ) dx 0 (3.14) =( e (( a 1) +1)+ +1 2 e e a ) (3.15) J L;p m +2 = Z 1 1 a (1 x 0 ) e (2 a x 0 ) dx 0 (3.16) =( e ( 1)+ e e a ) (3.17) J R;p m 1 = Z a 1 (1+ x 0 ) e ( x 0 ( 1 a ) dx 0 (3.18) =( e (( a 1) 1)+ e a ) (3.19) J R;p m = Z 1 a a (1 j x 0 j ) e ( x 0 ( a ) dx 0 (3.20) =((1 a ) +1+ e ( +1) 2 e a ) (3.21) J R;p m +1 = Z 1 1 a (1 x 0 ) e ( x 0 (1 a )) dx 0 (3.22) =( 1+ e a ) (3.23) Notethatjustoneevaluationofanexponentialfunctionisrequiredperparticle(namely e a ).Toaccountforarbitrary x ,wemakethesubstitution x = . Toobtainthetotallocaldeposit,wesimplysumtheparticlecontributionsontothegrid. Let N p bethetotalnumberofparticles.Thealgorithmforthelocaldepositstepisgiven by: Initialize J L k = J R k =0forall k . for p =1: N p do 68 Particle p locatedincell[ x m ;x m +1 ]. for j = m r : m + r do Compute J L;p j +1 , J R;p j Deposit J L j +1 = J L j +1 + J L;p j +1 , J R j = J R j + J R;p j endfor endfor Notethatthelocaldepositstepcosts O ( N p )operations. 3.4.1.2RecursiveSweepStepin1D Oncewehaveperformedthelocaldepositstep,wecompletetheevaluationoftheparticlein- tegralwitharecursivesweepstep.Supposewehave N gridpoints, x 1 ;:::;x N .Thealgorithm fortherecursivesweepstepisgivenby: Initialize I L 1 = I R N =0 for j =1: N 1 do I L j +1 = J L j +1 + e I L j I R N j = J R N j + e I R N j +1 endfor I = I L + I R Notethattherecursivesweepstepcosts O ( N )operations. 69 3.4.2FastConvolutionAlgorithmin2D Foraseparableparticleshape S ( x;y )= S x ( x ) S y ( y ),wehave I x [ I y [ S ]]( x;y )= I x [ S x ]( x ) I y [ S y ]( y )(3.24) = I L x [ S x ]( x )+ I R x [ S x ]( x ) I D y [ S y ]( y )+ I U y [ S y ]( y ) (3.25) = I L x [ S x ]( x ) I D y [ S y ]( y )+ I L x [ S x ]( x ) I U y [ S y ]( y )+(3.26) + I R x [ S x ]( x ) I D y [ S y ]( y )+ I R x [ S x ]( x ) I U y [ S y ]( y )(3.27) Thisissuggestiveofhowwewillbuildthe2Dalgorithm. 3.4.2.1LocalDepositStepin2D Considerparticle p centeredat( x p ;y p ) 2 [ x m ;x m +1 ] [ y n ;y n +1 ],withseparableparticle shape S p ( x;y )= S ( x x p x ) S ( y y p y )where S ( x )= 8 > > > < > > > : 1 j x jj x j < 1 ; 0 j x j 1 (3.28) ThesupportoftheparticleshapeisshowninFigure3.2. Inthelocaldepositstep,weformatensorproductonthegridassuggestedbytheabove decomposition.Notethatatotalof12localintegralsmustbeevaluatedforeachparticle, thensummedontothegridasatensorproduct. Initialize J LU j;k = J LD j;k = J RU j;k = J RD j;k =0forall j;k . for p =1: N p do Particle p locatedincell[ x m ;x m +1 ] [ y n ;y n +1 ]. 70 1 a x m 1 1 a x m 0 1 a x m +1 1 2 a x m +2 Figure3.1: S ( x )=1 j x j ; j x j < 1 ;S ( x )=0 ; j x j 1 x m 1 x m x m +1 x m +2 y n 1 y n y n +1 y n +2 ( x p ;y p ) Figure3.2:Thesupportofalinearparticleshape Sp ( x;y )in2D. 71 for j = m r : m + r , k = n r : n + r do Compute J L;p j +1 , J R;p j , J D;p k +1 , J U;p k Deposit: J LD j +1 ;k +1 = J LD j +1 ;k +1 + J L;p j +1 J D;p k +1 J LU j +1 ;k = J LU j +1 ;k + J L;p j +1 J U;p k J RD j;k +1 = J RD j;k +1 + J R;p j J D;p k +1 J RU j;k = J RU j;k + J R;p j J U;p k endfor endfor 3.4.2.2RecursiveSweepStepin2D Therecursivesweepstepissimilartothe1Dcase,andisgivenbelow. for k =1: N y do Initialize I LD 1 ;k = I LU 1 ;k = I RD N x ;k = I RU N x ;k =0. for j =1: N x 1 do I LD j +1 ;k = J LD j +1 ;k + e I LD j;k I LU j +1 ;k = J LU j +1 ;k + e I LU j;k I RD N x j;k = J RD N x j;k + e I RD N x j +1 ;k I RU N x j;k = J RU N x j;k + e I RU N x j +1 ;k endfor endfor J D = I LD + I RD J U = I LU + I RU for j =1: N x do 72 Initialize I D j; 1 = I U j;N y =0. for k =1: N y 1 do I D j;k +1 = J D j;k +1 + e I D j;k I U j;N y k = J U j;N y k + e I U j;N y k +1 endfor endfor I = I D + I U Both1Dand2Doverallalgorithmscost O ( N p + N )operations,where N isthetotal numberofgridpoints.SinceinatypicalPICsimulation, N p >>N ,thecostoftheoverall algorithmisdominatedbythelocaldepositstep. 3.5ParticleBoundaryConditions Indealingwithboundaries,twotypesofconsiderationsmustbemade.First,wemust determinewhattodointheintegrationofboundaryparticles,forwhichthesupportofthe shapefunctionextendsoutsideofthedomain.Thiswillbedependentuponthetypeof boundarycondition.Forperiodicboundaryconditions,wecansimplyextendtheparticle shapefunctionperiodicallyandproceedtointegrate.ForDirichletboundaryconditions,we extendtheintegrationdomaintoincludeghostpointsjustbeyondtheboundarytowhich theboundaryparticlesareweighted. Second,wemustensurethattheparticleconvolutionintegralisconsistentwiththe boundaryconditionsonthewavefunction.Thisiseasilyhandledthroughtheusualbound- arycorrectiontermsinone-dimension,andcanbeextendedtothedimensionally-splitmul- tidimensionalcase. 73 3.5.11DPeriodicBC Inthecaseof1Dperiodicboundaryconditions,weextendtheshapefunctionoftheboundary particleasinFigure3.3,andthelocaldepositstepiscarriedoutforthisextendedshape function.Therecursivesweepstepandtheboundarycorrectionsteparecarriedoutas usual. x N 1 x 1 x 2 x p x 3 x N 1 x N x 2 x 3 Figure3.3:Periodicextensionofshapefunctionforboundaryparticle. 3.5.22DPeriodicBC Periodicboundaryconditionsareeasilyimposedinhigherdimensionsbysimilarlyperi- odicallyextendingtheparticleshapefunctionsofboundaryparticles,peformingthelocal depositstepaccordinglyandtherecursivesweepstepasusual. 3.5.31DDirichletBC InimposinghomogeneousDirichletboundaryconditionswithboundaryparticles,wecon- sidertwopossibleapproaches.Theapproachsimplyextendstheintegrationdomainto includeappropriateghostcellsextendingjustpasttheboundarytoincludetheentiresupport oftheparticles.Aboundarycorrectiontermoftheform Ae x (attheleftboundary x =0) toimposetheboundarycondition.Consideringtheleftboundaryonly,weobtainatotal potentialoftheform x )= p ( x ) p (0) e x ,where p isthepotentialassociatedwith 74 theparticle.Thesecondapproachplacesappropriateimageparticlesontheoppositeside oftheboundary,whileappropriatelyextendingtheintegrationdomainwithghostcellsto includethesupportsoftheboundaryparticlesandtheirimageparticles.Consideringagain theleftboundaryonly,weobtainatotalpotentialoftheform x )= p ( x )+ img ( x ),where p isagainthepotentialassociatedwiththeparticleand img isthepotentialassociated withtheimageparticle. WecanverifythatbothapproachescanimposethehomogeneousDirichletboundary condition=0,butthepotentialandtheassociatedwiththeapproachesarenot identical.Duetothethatarisenearcornersorothergeometricsingularitiesof theboundary,itispreferabletoavoidtheuseofimageparticles.Wedescribetheabove approachesforlinearparticleshapes,buttheanalysiscanbeextendedtootherparticle shapes. 3.5.3.1AnalysisoftheFirstApproach Consideragrid x 1 =0, x j +1 = x j + x , j =1 ; 2 ;::: .Considerasingleboundaryparticle centeredinthecellat x p 2 [ x 1 ;x 2 ]=[0 ; x ].Theshapefunctionoftheparticleis S p ( x )= 8 > > > < > > > : 1 j x x p j = x j x x p j < x 0 j x x p j x (3.29) 75 Thepotentialassociatedwiththeparticleis p ( x )= I [ S p ]( x )(3.30) = Z x p x x p x e j x y j S p ( y ) dy (3.31) = Z x x p x e ( x y ) S p ( y ) dy +(3.32) + Z x p x x e ( y x ) S p ( y ) dy (3.33) = I L [ S p ]( x )+ I R [ S p ]( x )(3.34) Theassociatedwiththeparticleis E p ( x )= d dx p ( x )(3.35) = I L [ S p ]( x ) I R [ S p ]( x ) j x x p j < x (3.36) = I L [ S p ]( x ) x x p + x (3.37) ToimposeahomogeneousDirichletboundaryconditionat x =0withacorrectionterm Ae x ,wechoose A = p (0).Thus,thetotalpotentialintheapproachisgivenby x )= p ( x ) p (0) e x (3.38) 76 Thetotalfortheapproachisgivenby E ( x )= d dx x )(3.39) = I L [ S p ]( x ) I R [ S p ]( x )+ p (0) e x j x x p j < x (3.40) = I L [ S p ]( x )+ p (0) e x x x p + x (3.41) = E p ( x )+ p (0) e x (3.42) 3.5.3.2AnalysisoftheSecondApproach Considerasingleboundaryparticle,asbefore.Inthesecondapproach,weplaceanimage particlecenteredat x p .Theshapefunctionoftheimageparticleis S img ( x )= 8 > > > < > > > : 1+ j x + x p j = x j x + x p j < x 0 j x + x p j x (3.43) Thepotentialassociatedwiththeimageparticleis img ( x )= I [ S img ]( x )(3.44) = Z x p x x p x e j x y j S img ( y ) dy (3.45) Itiseasytoseethroughsymmetrythat img (0)= p (0).Ifweattempttoapplythe usualboundarycorrectionterm,thetotalpotential x )= p ( x )+ img ( x )+ Ae x must thensatisfy A =0.Sothetotalpotentialis x )= p ( x )+ img ( x )(3.46) 77 andtheboundarycondition=0. Theassociatedwiththeimageparticleis E img ( x )= d dx img ( x )(3.47) = I L [ S img ]( x ) I R [ S img ]( x ) j x + x p j < x (3.48) = I L [ S img ]( x ) x> x p + x (3.49) Thetotalisthen E ( x )= E p ( x )+ E img ( x ). 3.5.3.3ComparisonoftheApproaches Itisclearthattheinthepotentialbetweentheapproachesis x )= img ( x )+ p (0) e x (3.50) andtheceinthebetweentheapproachesis E ( x )= E img ( x ) p (0) e x (3.51) Wehaveanalyticformulasfortheseexpressionswecanuseforthecomparison.Wecan 78 computethat img ( x )= 8 > > > < > > > : 2 e ( x + x p ) x (cosh( x ) 1) x x p + x 2 1 ( x + x p ) j = x e ( x + x p ) x + e x x cosh( ( x + x p )) 0 x< x p + x (3.52) Ontheotherhand, p (0) e x =2 " 1 x p = x e x p x + e x x cosh( x p ) # e x (3.53) Fromtheseformulas,wecanshowthat x )=2 e x 1 x p x + 1 x sinh( x x p )) (3.54) = O (( x ) 2 )(3.55) andthat E x )= x )(3.56) = O ( 3 x 2 )(3.57) 79 3.5.42DDirichletBC Ifweusetheapproachdescribedabove,wewillextendtheintegrationdomainwith ghostcellsbeyondtheboundarytoincludethesupportofallparticles.Thentheboundary conditionwillbeimposedwiththeusualboundarycorrectionterm. 3.6NumericalResults 3.6.1ElectrostaticTestProblems Weconsiderthreestandardperiodicelectrostatictestproblemsin1Dand2D,then a1Dboundedplasmaproblem,thesimulationofsheathformation.Inthethreetest problems,electronsareloadedfromaperturbedinitialdistributionoftheform f e ( x;v;t =0)= f e ( v ) 1+ sin 2 ˇx L x (3.58) where L x isthelengthofthedomain, istheamplitudeofperturbation,and f e ( v )isthe initialvelocitydistribution.Inthe2Dcase,simulationsaretakentobeuniforminthe y - direction.Wenormalizequantitiesaccordingtothenondimensionalizationpresentedinthe Section1.3.1.Inparticular,wenormalizetimequantitiestotheinverseplasmafrequency ! 1 p .Wewillconsideraperiodicdomainwithauniformneutralizingbackgroundcharge, andfurtherwesetthespeedoflight c =100.Intheseproblems,weseegoodperformance evenatlargeCFLs,sincethephysicsisdominatedbythelowfrequencyspatialmodes. 80 3.6.1.1ColdPlasmaLangmuirWave WeconsideracoldplasmaLangmuirwave[4]with f e ( v )= ( v ).Electronsareperturbed awayfromauniformlydistributed,motionlessstateagainstastatic,uniformneutralizing backgroundchargedistribution.Theresultingseparationofchargeproducescoldplasma oscillation.Inthe1Dsimulation,weset L x =2 ˇ and =0 : 1.Weusea100cellgridand take t =0 : 1,andweuse N p =3600particles.Inthe2Dsimulation,weset L x = L y =2 ˇ andtheperturbationstrength =0 : 1.Weusea100 100gridandagaintake t =0 : 1, andweuse N p =360000particles.Inboth1Dand2Dcases,theCFLnumberusedis c t= x ˇ 159,muchlargerthanwhatwouldbeallowedbyanexplicitmethod.The oscillationinthepotentialenergyisplottedandcomparedtothepredictionoflineartheory inFigure3.6;weseethattheplasmafrequencyisaccuratelyreproduced. 3.6.1.2TwoStreamInstability Weconsiderthetwostreaminstabilitywith f e ( v )= ( v v beam )+ ( v + v beam ).Two counterstreamingbeamsofelectronsareperturbedawayfromauniformlydistributedstate againstastatic,uniformneutralizingbackgroundchargedistribution.Thebeamsinteract and\rollup"inphasespace,causingsomeoftheparticles'kineticenergytobetransformed intopotentialenergystoredintheelectricAccordingtothedispersionrelationforthe twostreaminstabilityfromlineartheory[4],wehave ! 4 2 ! 2 ( ! 2 p + k 2 v 2 beam )+ k 2 v 2 beam ( k 2 v 2 beam 2 ! 2 p )=0(3.59) whichgivesthegreatestgrowthrate, ˇ 0 : 3535,for k ˇ 3 : 06.Wethereforescalethe domaintothisvalueof k ,andtake L x =2 ˇ= 3 : 06.Wetakethebeamvelocity v beam =1 81 x 0 x 1 x 2 x p x 3 Figure3.4:Particleshapefunction S p ( x )ofaboundaryparticlewithghostcell[ x 0 ;x 1 ]= [ x; 0]. x 1 x 0 x p x 1 x 2 x p x 3 Figure3.5:Particleshapefunctions S p ( x )(solid)and S img ( x )(dotted)ofaboundaryparticle anditsimageparticle,withghostcells[ x 1 ;x 0 ]=[ x; x ]and[ x 0 ;x 1 ]=[ x; 0]. 10 20 30 40 50 60 0 0 : 2 0 : 4 0 : 6 0 : 8 1 Time ( 1 =! p ) PotentialEnergy Figure3.6:Potentialenergyincoldplasmaoscillation.Greenisthe1Dnumericalresult, blueisthe2Dnumericalresult,andredisthepredictionoflineartheory.Weseetheplasma frequencyisaccuratelyreproducedinoursimulations. 82 andtheperturbationstrength =0 : 0005.Inour1Dsimulation,weusea100cellgrid with t =0 : 1,andweuse N p =1000particles.Inour2Dsimulation,weusea100 100 gridwith t =0 : 1,andweuse N p =1000000particles.ThisresultsinaCFLnumberof c t= x ˇ 68.Werunthesimulationsfor1000timesteps.Thegrowthofthe k =3 : 06 modeoftheelectricisshowninFigure3.7forthe1Dand2Dcases,andagreeswith theratefromlineartheory.Inthenonlinearsaturationstage,weseeaslightdiscrepancy betweenthe1Dand2Dresults,probablyduetotheaccumulationofnumericalerror.We alsoshowselectedphasespaceplotsinFigure3.8,whereweseetheexpected\rollingup" ofthetwobeams.Resolutionislimitedbythenumberofparticles. 0 10 20 30 40 50 60 70 80 90 100 10 5 10 3 10 1 10 1 10 3 Time ( 1 =! p ) ElectricField-LowestModeCo 1Dresult 2Dresult(y=0slice) lineartheory Figure3.7:Growthofthemodewithmaximumgrowthrateinthetwostreaminstability, correspondingto k =3 : 06.Greenisthe1Dnumericalresult,blueisthe2Dnumericalresult (measuredalongthecentral y =0slice),andredisthepredictionoflineartheory.Wesee thecorrectgrowthrateisreproducedinoursimulations. 83 (a) (b) (c) (d) (e) (f) (g) (h) Figure3.8: Weseeselectedparticlephasespaceplotsforthetwostreaminstabilityproblem.The leftcolumnisfromthe1Dsimulation,whiletherightcolumnisfromthe2Dsimulation,following asliceofparticlesinitializedalongtheline y = L y = 2. 84 3.6.1.3LandauDamping WeconsiderLandaudampingofLangmuirwavesinawarmplasma,with f e ( v )takento beMaxwellian.Warmelectrons,followingaMaxwellianvelocitydistribution,areperturbed awayfromaunifomdistributionagainstastatic,uniformneutralizingbackgroundcharge distribution.Potentialenergyfromtheelectricistransformedintokineticenergyof particles.Thedispersionrelationfromlineartheoryinthiscasegivesadecayrateof ˇ 0 : 154forthe k =0 : 5mode[4].Wetake L x =4 ˇ ,electronthermalvelocity v therm =1 andperturbationstrength =0 : 1.Inour1Dsimulation,weusea100cellgridandtake t =0 : 1,andweuse N p =1000000particles.Inour2Dsimulation,weusea100 100 gridandtake t =0 : 1,andweuse N p =9000000particles.Werunthesimulationsfor300 timesteps.Thedecayofthe k =0 : 5modeoftheelectricinthe1Dand2Dsimulations isshowninFigure3.9,andagreeswiththeratefromlineartheory.Asinthetwostream instabilityexample,thereisadiscrepancybetweenthe1Dand2Dresultsatlatertimes, againlikelyduetotheaccumulationofnumericalerrors. 3.6.1.4SheathFormationinaBounded1DPlasma Wepresentthesimulationofsheathformationinabounded1Dplasma,followingthemodel describedin[46].Incontrasttothepreviousproblems,thissimulationincorporatesboth mobileelectronsandions.ElectronsandionsareinitializedfromMaxwelliandistributions anduniformlyspatiallydistributedinaboundeddomain.Theleftboundaryisasymmetry plane,andsoweimposeNeumannboundaryconditionsonthepotential,andboundary conditionsonparticles,asin[46].Therightboundaryisaconductorthatcollectscharged particles.Whenparticleshittherightboundary,theyareremovedfromthesimulation. Sinceelectronshaveahigheraveragevelocitythanions,theyhaveagreateronthe 85 collectorandbecomedepletedneartherightboundary,wherethencebetweenion andelectrondensitiesleadstoacollectorsheathregion,wherethepotentialchangesfrom theinteriorvaluetothewallvalue,whichhastheofrepellingelectronsawayfrom therightwall.Electronsandionsarereplenishedbyaparticlesourceregionneartheleft boundary,whereelectronsandionsareinjecteduniformlyintheregionfromaMaxwellian distributionataratepertimestep.Wetake m i =m e =100, T src;i =T src;e =1,and weset v th;e =1,andset L x =20(inDebyelengths).Weusea100cellgridandtake t =0 : 1,whichgivesaCFLnumberof50.Werunoursimulationfor8000timesteps, upto3 : 6thermal-iontransittimes.InFigure3.10weseetheresultofthesimulation.In Figure3.10a,weseetheofthepotential,whichhastherightqualitativefeatures, includingacollectorsheathregionthatisseveralDebyelengthswide.InFigure3.10b,we seethenetelectronandioncounts,alongwiththeinjectionrate.Thebetween theelectronandioncountstheencebetweenelectronandiondensitiesinthe collectorsheathregion. 3.6.2ElectromagneticTestProblems 3.6.2.1BennettPinchProblem WepresenttheapplicationofourPICmethodtotheBennettpinch[3],anrelated tothemagnetictofabeamofchargedparticles.Abeamofchargedparticles inducesasolenoidalmagneticaroundthebeam.Particlesneartheedgeofthebeam moveorthogonallytotheselinesatthebeamdriftvelocity,causingtheparticlestobe acceleratedtowardsthecenterofthebeam,inparticlesinthebeam.An appropriatechoiceofparametersleadstoastationarysteadystate,uniformalongtheaxisof 86 thebeam.Awell-knownmagnetohydrodynamic(MHD)modelofastationarysteadystate givesexplicitformulasforthebeamdensityandthemagneticeld[6],andprovidesabasis forthevalidationofournumericalmethod.Moreover,itisasteptowardapplyingour methodtomorephysicallyinterestingbeaminstabilityproblemsinthreedimensions. OurPICsimulationoftheBennettpinchistwo-dimensionalinphysicalspace,andthree- dimensionalinvelocityspace.Theparticlebeamisconsidereduniformalongtheaxisofthe beam,whichwetaketobethe z -direction,whichreducesthephysicaldimensionstotwo. Electronsdriftinthe z -directionwithauniformaveragebeamdriftvelocity v b ,andthis motioninducesamagneticwithonly x -and y -components.Astationaryion backgrounddistributionenforcesquasineutralityinthebeam,withanyseparationofcharge producinganelectricwithonly x -and y -componentsactingasarestoringforce.The electronsareassumedtofollowaMaxwelliandistributionwiththermalvelocity v th .We take v b =v th =100and c=v th =1000,where c isthespeedoflight.Theionsareconsidered cold( T i =0). Sincethebeamdriftvelocityistakentobemuchlargerthanthe(transverse)thermal velocity,andfurther,thetransversevelocitiesfollowaMaxwelliandistributionandsoshould notgenerateanynetcurrents,weneglectthe x -and y -componentsofthecurrentdensity(and soalsoof A ).Inthetruesolution, @ @z =0and @A z @t =0,soweneglect E z = @ @z @A z @t . Hence,weactuallyonlysolvetwowaveequation,onefor A z ,obtainingonlytransverse magneticcomponents, B x = @A z @y and B y = @A z @x ,andoneforobtainingonly transverseelectriccomponents E x = @ @x and E y = @ @y .Thus,thePoissonequation bythescalarpotentialis = 0 .Wediscretizeourdomainwithastaggered grid,onecellofwhichisshowninFigure3.11. Thescalarpotentialiscalculatedfromthestandard5-pointLaplacian, 87 andtheequation i +1 ;j + i;j +1 + i 1 ;j + i;j 1 i;j x 2 = ˆ i;j 0 : (3.60) Theelectricandmagneticsarecalculatedonthestaggeredgridby as E i +1 = 2 ;j x = i +1 ;j i;j x (3.61) E i;j +1 = 2 y = i;j +1 i;j x (3.62) B i;j +1 = 2 x = A i;j +1 z A i;j z x (3.63) B i +1 = 2 ;j y = A i +1 ;j z A i;j z x : (3.64) TheelectricdthenthefollowingdiscreteanalogueofGauss'law: [ r E ] i;j = E i +1 = 2 ;j x E i 1 = 2 ;j x x + E i;j +1 = 2 y E i;j 1 = 2 y x (3.65) = 1 x i +1 ;j i;j x ! i;j i 1 ;j x ! +(3.66) + i;j +1 i;j x ! i;j i;j 1 x !! (3.67) = i +1 ;j + i;j +1 + i 1 ;j + i;j 1 i;j x 2 (3.68) = ˆ i;j 0 : (3.69) Themagneticthefollowingdiscreteanalogueofthedivergencefreecondi- 88 tion: [ r B ] i +1 = 2 ;j +1 = 2 = B i +1 ;j +1 = 2 x B i;j +1 = 2 x x + B i +1 = 2 ;j +1 y B i +1 = 2 ;j y x (3.70) = 1 x A i +1 ;j +1 z A i +1 ;j z x ! A i;j +1 z A i;j z x ! +(3.71) + A i +1 ;j +1 z A i;j +1 z x ! A i +1 ;j z A i;j z x !! (3.72) =0 : (3.73) Allcomputationalboundariesinthisproblemarewboundaries.Inordertosupply thePoissonsolverwithsuitableboundaryvalues,thewavesolverisapplied withwboundariesconditionstoevolvethewavepotential W alongsidethePoisson potentialTheboundaryvaluesfrom W arethensuppliedtothePoissonsolvertouse incalculatingOncethewavesolverreachessteadystate, W andonlyby0.1% relativeerror,however,thewavepotentialgivesadiscreteelectricwithdivergenceerror ontheorderof10 3 ,whilethePoissonpotentialgivesadiscreteelectricwithdivergence errorontheorderofmachineepsilon10 16 . Likeintheothertestproblems,weusetheeversionofwavesolver.Particle velocitiesinallthreedirectionsareupdatedwiththenonrelativisticBorispush[4].Particles areinitializedaccordingtotheMHDsteadystate(accordingtothetheoreticalspatialdensity andthecorrespondingMaxwelliandistributioninvelocityspace)andheldwhile thesolverissteppedtoanapproximatesteadystate,afterwhichtheparticlepushis turnedon.Thesimulationisruntoatimeof R b = (2 v th )(plusstartuptime),where R b isanebeamradiusand v th isthethermalvelocity,atwhichtimetherewould 89 besubstantialspreadingofthebeaminabsenceofthecotWechoose R b suchthat99%oftheparticlesinthetheoreticalbeamarewithinthisradius.Inloading particles,thebeamiscutatradius R b (noparticlesareloadedoutsideofthisradius). Thecomputationaldomainistakentobea4 R b 4 R b squarecenteredonthebeamaxis.The computationaldomainistruncatedwithwboundaryconditions.Particlesexitingthe boundaryofthecomputationaldomainarereinjectedintothebeamtomaintainconstant totalcurrent.However,sincemostparticlesshouldbewithinthebeam,such boundarycrossingsshouldberare. NumericalresultsfortheBennettpincharegiveninFigure3.12.Inordertoresolve largegradientsnearthecenterofthebeam,weusea500-cellby500-cellgridandaCFL numberof c t= x =3(exceptinFigure3.12casnoted)and500,000electronparticles. Finalnumericalsolutionsareshownatthetime,after334startuptimestepsand 20,834PICtimesteps.(Thetimeisapprox.22plasmaperiods,andthediameterof thebeamisapprox.280Debyelengths.)InFigure3.12a,weseegoodagreementbetween thenumericalelectrondensityandMHDtheory.Theinsetzoomedportionshowsaslight discrepencyatthepeakofthebeam,duetostatisticaltuationcausedbythe numberofparticles.InFigure3.12b,weseethetimehistoriesofthepotentialenergy, calculatedas P j x y 2 h 1 R ( B 2 x;j + B 2 y;j )+ " R ( E 2 x;j + E 2 y;j ) i wherethesumisovergrid points j (with " R and R asinSection1.3.2),andthekineticenergy,calculatedas P i 1 2 m i ( v 2 x;i + v 2 y;i +( v x;i v b ) 2 )wherethesumisoverelectronparticles i .Weseegood energyconservation,despitetheslightdissipationoftheescheme.Theinitialspike inthepotentialenergyistheresultoftransientwaves,arisingduetothebeamturningon, andwingoutofthedomainasthesolutionissteppedtoasteadystate.InFigure3.12c, weseetheresultoftin t ,keeping x showingapoftheazimuthal 90 magnetic B alongthecentral y =0sliceforCFLnumbersof3,10and20,alongwith MHDtheory.Weobserveapproximatesecond-orderconvergencein t ,asexpected(amore robustconvergencestudyisconfoundedbytheslowconvergenceinparticlenumberinPIC methods).Outsideofthebeamradius R b =1,thereiserrorassociatedwiththe radiusofthebeam(thetheoreticalbeamdensitydecaysonlyalgebraically).InFigure3.12d, weseethenumericalerroroftheazimuthalmagnetic B (withaCFLnumberofof 3)normalizedbythepeakvalueofthemagneticandweseethatthereisageometric patterntothenumericalerror,characteristicofthedimensionally-splitmethod.Inaddition tothissplittingerror,thetotalerroriscontributedtobyerrorsassociatedwiththespatial quadratureandtheusedtocalculatethemagnetic(likelycontributing tothelargeerroratthecenterofthebeamduetolargegradientsthere)andwiththe beamradiusandthewboundarycondition(contributingmoststronglynear theboundaryofthecomputationaldomain).Theseresultsshowthatourmethodcanindeed simulateabasicelectromagneticplasmaphenomenonwithaCFLnumberlargerthanwhat isallowedbytypicalexplicitschemes.TheCFLnumberusedinthisproblemislimited bytheaccuracyofthesecond-orderwavesolver.Ahigher-orderwavesolver,suchasthose in[11],wouldallowforalargerusabletimestepsize,andwillbethesubjectoffurther investigation. 3.6.2.2MardahlBeamProblem Weapplyourmethodtothebeamproblemproposedin[37]asadiagnosticforthe ofdivergenceerror.Particlesareinjectedintothedomain,whichisaboxwithPECwalls, fromtheleftwall,travelacrossthedomainandareremovedfromthesimulationastheyhit therightwall.Parametersarechosensuchthatthebeamshouldpassthroughthedomain 91 unperturbed. IntheMardahlbeamproblem,wehavecurrentsintheplaneofsimulationonly,andso weretainthe x -and y -componentsofthevectorpotential, A x and A y ,alongwiththescalar potentialinthemodel.Weretaintheelectriccomponents E x = @ @x @A x @t and E y = @ @y @A y @t andthemagneticcomponent B z = @A y @x @A x @y .ThePoissonequation bythescalarpotentialis = 0 + @ @t @A x @x + @A y @y (3.74) Wediscretizeourdomainwithastaggeredgrid,onecellofwhichisshowninFigure3.13. Denotingby D t a(linear)discretizationofthetimederivativeoperator @=@t ,thescalarpotentialtheequation i +1 ;j + i;j +1 + i 1 ;j + i;j 1 i;j x 2 = ˆ i;j 0 +(3.75) D t A i +1 = 2 ;j x A i 1 = 2 ;j x x + A i;j +1 = 2 y A i;j 1 = 2 y x ! : (3.76) Theelectricandmagneticsarecalculatedonthestaggeredgridby 92 as E i +1 = 2 ;j x = i +1 ;j i;j x D t ( A i +1 = 2 ;j x )(3.77) E i;j +1 = 2 y = i;j +1 i;j y D t ( A i;j +1 = 2 y )(3.78) B i +1 = 2 ;j +1 = 2 z = A i +1 ;j y A i;j y x A i;j +1 x A i;j x y (3.79) TheelectricdthenthefollowingdiscreteanalogueofGauss'law: [ r E ] i;j = E i +1 = 2 ;j x E i 1 = 2 ;j x x + E i;j +1 = 2 y E i;j 1 = 2 y y (3.80) = 1 x i +1 ;j i;j x D t ( A i +1 = 2 ;j x ) ! i;j i 1 ;j x D t ( A i 1 = 2 ;j x ) ! + (3.81) + i;j +1 i;j x D t ( A i;j +1 = 2 y ) ! i;j i;j 1 x D t ( A i;j 1 = 2 y ) !! (3.82) = i +1 ;j + i;j +1 + i 1 ;j + i;j 1 i;j x 2 (3.83) D t A i +1 = 2 ;j x A i 1 = 2 ;j x x + A i;j +1 = 2 y A i;j 1 = 2 x x ! (3.84) = ˆ i;j 0 : (3.85) wherewehaveusedthelinearityof D t . NumericalresultsforthisproblemaregiveninFigure3.14,usinga64 64gridand aCFLof1.Weseetheexpecteddistortionofthebeaminthecasewhenthedivergence errorisnotcontrolledthroughtheellipticcorrection,whereasthebeampassesthroughthe 93 domainunperturbed,asexpected,whenusingthePoisson-basedpotential. 94 0 5 10 15 20 25 30 10 5 10 4 10 3 10 2 10 1 10 0 Time ( 1 =! p ) ElectricField-LowestModeCo 1Dresult 2Dresult(y=0slice) lineartheory Figure3.9:Landaudampingofthelowestmode,correspondingto k =0 : 5.Greenisthe1D numericalresult,blueisthe2Dnumericalresult(measuredalongthecentral y =0slice), andredisthepredictionoflineartheory.Weseethatthecorrectdecayrateisreproduced inoursimulations 95 0 2 4 6 8 10 12 14 16 18 20 0 : 8 0 : 6 0 : 4 0 : 2 0 0 : 2 SourceRegion Particle/NeumannBC(Sym.Plane) CollectorSheath ParticleCollector Position ( D ) ScalarPotential (a) 0 0 : 5 1 1 : 5 2 2 : 5 3 3 : 5 4 0 : 5 1 1 : 5 2 2 : 5 3 10 4 Time(thermal-iontransittimes) NumberofParticles (b) Figure3.10: In3.10a,weseethescalarpotentialleat t =3 : 6thermal-iontransittimes.In 3.10b,weseethesimulationelectronandioncount,theredandbluecurvesrespectively,along withtheinjectionrate,theblackdashedline. 96 i;j , A i;j z ˆ i;j , J i;j z E i +1 = 2 ;j x B i +1 = 2 ;j y E i 1 = 2 ;j x B i 1 = 2 ;j y E i;j +1 = 2 y , B i;j +1 = 2 x E i;j 1 = 2 y , B i;j 1 = 2 x Figure3.11:ThestaggeredgridusedfortheBennettpinchproblem. 97 (a) (b) (c) (d) Figure3.12: Theshows3.12aelectrondensity,3.12bandpotentialenergies,3.12cmagnetic atvariousCFLnumbers,and3.12dtherelativeerrorintheazimuthalmagnetic B (normalizedtothemaximumvalueofthemagneticResultsarewithaCFLnumberof3, exceptasnotedin3.12c.Positionunitsareintermsofebeamradius R b . 98 B i +1 = 2 ;j +1 = 2 z J i +1 ;j +1 = 2 y A i +1 ;j +1 = 2 y E i +1 ;j +1 = 2 y J i;j +1 = 2 y A i;j +1 = 2 y E i;j +1 = 2 y J i +1 = 2 ;j +1 x , A i +1 = 2 ;j +1 x , E i +1 = 2 ;j +1 x J i +1 = 2 ;j x , A i +1 = 2 ;j x , E i +1 = 2 ;j x ˆ i;j , i;j ˆ i +1 ;j , i +1 ;j ˆ i;j +1 , i;j +1 ˆ i +1 ;j +1 , i +1 ;j +1 Figure3.13:ThestaggeredgridusedfortheMardahlbeamproblem. 99 (a) (b) (c) (d) Figure3.14: Theshowsthedivergenceerrorintheelectricandthebeamdistri- butioncalculatedfromawaveequationpotential3.14a,3.14candinthePoissonequationpotential 3.14b,3.14d. 100 Chapter4 Conclusion Inthiswork,wehavedescribedaPICmethodthatusesanunconditionallystablewave equationsolvertoeliminatetheCFLrestrictionontheratioofthetimestepsizetothe spatialstepsize,typicalofexplicitmethods,whileretainingcomputationalcostandcode complexitycomparabletosuchexplicitmethods.Ournumericalresultsshowthatwecan applyourmethodtoproblemsofplasmaphysicsusingatimestepsizelargerthanwhat wouldbeallowedbyatypicalexplicitsolver.Wehaveseenthattheusabletimestep sizecanbelimitedbythenumericalaccuracyofthemethodwhentherearelargegradients (high-frequencycontent)inthesolution.Infuturework,wewillinvestigatetheuseofhigher- ordermethods,suchasgivenin[11],inourPICmethodinordertoincreasethemaximum usabletimestepsize,andwewillmakeuseoftheimplicitwavesolversabilitytohandle complexboundarygeometrieswithouttheuseofastaircasingapproximation.Afurther courseofactionwillbetoimplementaboundaryintegraltreecode(BIT)solutiontosolve themoHelmholtzequationsinthesemi-discreteschemes,suchasin[35],ratherthan usedimensionalsplitting. 101 BIBLIOGRAPHY 102 BIBLIOGRAPHY [1] TDArberandRGLVann, Acriticalcomparisonofeulerian-grid-basedvlasovsolvers , Journalofcomputationalphysics 180 (2002),no.1,339{357. [2] KiyoshiAsanoandSeijiUkai, Onthevlasov-poissonlimitofthevlasov-maxwellequation , StudiesinMathematicsandItsApplications 18 (1986),369{383. [3] WillardHBennett, Magneticallyself-focussingstreams ,PhysicalReview 45 (1934), no.12,890. [4] CharlesKennedyBirdsallandAllanBruceLangdon, Plasmaphysicsviacomputersim- ulaition ,CRCPress,2005. [5] GeorgeDavid Generalmeanvalueandremaindertheoremswithapplicationsto mechanicalentiationandquadrature ,TransactionsoftheAmericanMathematical Society 7 (1906),no.1,107{136. [6] JoseABittencourt, Fundamentalsofplasmaphysics ,Springer,2004. [7] JUBrackbillandDWForslund, Animplicitmethodforelectromagneticplasmasimu- lationintwodimensions ,JournalofComputationalPhysics 46 (1982),no.2,271{308. [8] OscarPBrunoandMarkLyon, High-orderunconditionallystablefc-adsolversfor generalsmoothdomainsi.basicelements ,JournalofComputationalPhysics 229 (2010), no.6,2009{2033. [9] M.Causley,Y.Gu˘clu,E.Wolf,andA.Christlieb, AFast,Unconditionallystablesolver forthewaveequationbasedontheMethodofLinesTranspose , (inpreparation) . [10] MatthewCausley,AndrewChristlieb,BenjaminOng,andLeeVanGroningen, Method oflinestranspose:Animplicitsolutiontothewaveequation ,MathematicsofCompu- tation(2014). [11] MatthewFCausleyandAndrewJChristlieb, Higherordera-stableschemesforthe waveequationusingasuccessiveconvolutionapproach ,SIAMJournalonNumerical Analysis 52 (2014),no.1,220{235. 103 [12] RomanChapko,RainerKress,etal., Rothesmethodfortheheatequationandboundary integralequations ,JournalofIntegralEquationsandApplications 9 (1997),no.1,47{69. [13] GuangyeChen,Luison,andDanielCBarnes, Anenergy-andcharge-conserving, implicit,electrostaticparticle-in-cellalgorithm ,JournalofComputationalPhysics 230 (2011),no.18,7018{7036. [14] Chio-ZongChengandGeorgKnorr, Theintegrationofthevlasovequationinc rationspace ,JournalofComputationalPhysics 22 (1976),no.3,330{351. [15] YingdaCheng,AndrewJChristlieb,andXinghuiZhong, Energy-conservingdiscontinu- ousgalerkinmethodsforthevlasov{maxwellsystem ,JournalofComputationalPhysics 279 (2014),145{173. [16] YingdaCheng,IreneMGamba,FengyanLi,andPhilipJMorrison, Discontinuous galerkinmethodsforthevlasov{maxwellequations ,SIAMJournalonNumericalAnalysis 52 (2014),no.2,1017{1049. [17] AndrewJChristlieb,RobertKrasny,JohnPVerboncoeur,JeroldWandIainD Boyd, Grid-freeplasmasimulationtechniques ,PlasmaScience,IEEETransactionson 34 (2006),no.2,149{165. [18] BruceICohen,ABruceLangdon,andAlexFriedman, Implicittimeintegrationfor plasmasimulation ,JournalofComputationalPhysics 46 (1982),no.1,15{38. [19] GaryCohen, Higher-ordernumericalmethodsfortransientwaveequations ,Springer, 2002. [20] SupriyoDeyandRajMittra, Alocallyconformalencetime-domain(fdtd) algorithmformodelingthree-dimensionalperfectlyconductingobjects ,Microwaveand GuidedWaveLetters,IEEE 7 (1997),no.9,273{275. [21] LawrenceC.Evans, Partialentialequations ,2ndeded.,Graduatestudiesinmath- ematics,vol.v.19,AmericanMathematicalSociety,Providence,R.I.,2010. [22] FrancisFilbetandEricSonnendrucker, Comparisonofeulerianvlasovsolvers ,Com- puterPhysicsCommunications 150 (2003),no.3,247{266. [23] MatthewRGibbonsandDennisWHewett, Thedarwindirectimplicitparticle-in- cell(dadipic)methodforsimulationoflowfrequencyplasmaphenomena ,Journalof ComputationalPhysics 120 (1995),no.2,231{247. 104 [24] DennisWHewettandABruceLangdon, Electromagneticdirectimplicitplasmasimu- lation ,JournalofComputationalPhysics 72 (1987),no.1,121{155. [25] RogerWHockneyandJamesWEastwood, Computersimulationusingparticles ,CRC Press,1988. [26] GBJacobsandJanSHesthaven, Implicit{explicittimeintegrationofahigh-order particle-in-cellmethodwithhyperbolicdivergencecleaning ,ComputerPhysicsCommu- nications 180 (2009),no.10,1760{1767. [27] GBJacobsandJSHesthaven, High-ordernodaldiscontinuousgalerkinparticle-in-cell methodonunstructuredgrids ,JournalofComputationalPhysics 214 (2006),no.1, 96{121. [28] Heinz-OttoKreiss,NAndersPetersson,andJacobom, enceapproximations oftheneumannproblemforthesecondorderwaveequation ,SIAMJournalonNumerical Analysis 42 (2004),no.3,1292{1323. [29] MaryCatherineAKropinskiandBryanDQuaife, Fastintegralequationmethodsfor rothe'smethodappliedtotheisotropicheatequation ,Computers&Mathematicswith Applications 61 (2011),no.9,2436{2446. [30] , Fastintegralequationmethodsforthemodhelmholtzequation ,Journalof ComputationalPhysics 230 (2011),no.2,425{434. [31] LevDavidovichLandauandEvgeniiMikhailovichLifshits, Theclassicaltheoryof , vol.2,Butterworth-Heinemann,1975. [32] ABruceLangdon, Onenforcinggauss'lawinelectromagneticparticle-in-cellcodes , ComputerPhysicsCommunications 70 (1992),no.3,447{450. [33] ABruceLangdon,BruceICohen,andAlexFriedman, Directimplicitlargetime-step particlesimulationofplasmas ,JournalofComputationalPhysics 51 (1983),no.1, 107{138. [34] JongwooLeeandBengtFornberg, Someunconditionallystabletimesteppingmethods forthe3dmaxwell'sequations ,JournalofComputationalandAppliedMathematics 166 (2004),no.2,497{523. [35] PeijunLi,HansJohnston,andRobertKrasny, Acartesiantreecodeforscreenedcoulomb interactions ,JournalofComputationalPhysics 228 (2009),no.10,3858{3868. 105 [36] MarkLyonandOscarPBruno, High-orderunconditionallystablefc-adsolversforgen- eralsmoothdomainsii.elliptic,parabolicandhyperbolicpdes;theoreticalconsiderations , JournalofComputationalPhysics 229 (2010),no.9,3358{3381. [37] PJMardahlandJPVerboncoeur, Chargeconservationinelectromagneticpiccodes; spectralcomparisonofboris/dadiandlangdon-mardermethods ,Computerphysicscom- munications 106 (1997),no.3,219{229. [38] StefanoMarkidisandGiovanniLapenta, Theenergyconservingparticle-in-cellmethod , JournalofComputationalPhysics 230 (2011),no.18,7037{7052. [39] MartinMasekandPaulGibbon, Mesh-freemagnetoinductiveplasmamodel ,Plasma Science,IEEETransactionson 38 (2010),no.9,2377{2382. [40] RLMorseandCWNielson, Numericalsimulationoftheweibelinstabilityinoneand twodimensions ,PhysicsofFluids(1958-1988) 14 (1971),no.4,830{840. [41] C-DMunz,PascalOmnes,RudolfSchneider,EricSonnendrucker,andUrsulaVoss, Di- vergencecorrectiontechniquesformaxwellsolversbasedonahyperbolicmodel ,Journal ofComputationalPhysics 161 (2000),no.2,484{511. [42] TakefumiNamiki, 3-dadi-fdtdmethod-unconditionallystabletime-domainalgorithm forsolvingfullvectormaxwell'sequations ,MicrowaveTheoryandTechniques,IEEE Transactionson 48 (2000),no.10,1743{1748. [43] DonaldWPeacemanandHenryHRachford,Jr, Thenumericalsolutionofparabolic andellipticentialequations ,JournaloftheSocietyforIndustrial&AppliedMath- ematics 3 (1955),no.1,28{41. [44] SVPetropavlovskyandSVTsynkov, Quasi-lacunaeofmaxwell'sequations ,SIAMJour- nalonAppliedMathematics 71 (2011),no.4,1109{1122. [45] EPohn,MagdiShoucri,andGKamelander, Eulerianvlasovcodes ,Computerphysics communications 166 (2005),no.2,81{93. [46] RJProcassini,CKBirdsall,andECMorse, Afullykinetic,self-consistentparticlesim- ulationmodelofthecollisionlessplasma{sheathregion ,PhysicsofFluidsB:Plasma Physics(1989-1993) 2 (1990),no.12,3191{3205. [47] Jing-MeiQiuandAndrewChristlieb, Aconservativehighordersemi-lagrangianweno methodforthevlasovequation ,JournalofComputationalPhysics 229 (2010),no.4, 1130{1149. 106 [48] JamesARossmanithandDavidCSeal, Apositivity-preservinghigh-ordersemi- lagrangiandiscontinuousgalerkinschemeforthevlasov{poissonequations ,Journalof ComputationalPhysics 230 (2011),no.16,6203{6232. [49] ErichRothe, Zweidimensionaleparabolischerandwertaufgabenalsgrenzfalleindimen- sionalerrandwertaufgaben ,MathematischeAnnalen 102 (1930),no.1,650{670. [50] JackSc Theclassicallimitoftherelativisticvlasov-maxwellsystem ,Communi- cationsinmathematicalphysics 104 (1986),no.3,403{421. [51] DavidNSmithe,JohnRCary,andJohanACarlsson, Divergencepreservationintheadi algorithmsforelectromagnetics ,JournalofComputationalPhysics 228 (2009),no.19, 7289{7299. [52] EricSonnendrucker,JeanRoche,PierreBertrand,andAlainGhizzo, Thesemi- lagrangianmethodforthenumericalresolutionofthevlasovequation ,Journalofcom- putationalphysics 149 (1999),no.2,201{220. [53] AlenTveandSusanCHagness, Computationalelectrodynamics:thefdtdmethod , ArtechHouseBoston,London(2000). [54] JohnPVerboncoeur, Aliasingofelectromagneticinstairstepboundaries ,Com- puterphysicscommunications 164 (2004),no.1,344{352. [55] , Particlesimulationofplasmas:reviewandadvances ,PlasmaPhysicsandCon- trolledFusion 47 (2005),no.5A,A231. [56] KaneSYee, Numericalsolutionofinitialboundaryvalueproblemsinvolvingmaxwell's equations ,IEEETrans.AntennasPropag 14 (1966),no.3,302{307. [57] FenghuaZhen,ZhizhangChen,andJiazongZhang, Towardthedevelopmentofathree- dimensionalunconditionallystableencetime-domainmethod ,Microwave TheoryandTechniques,IEEETransactionson 48 (2000),no.9,1550{1558. [58] FenghuaZheng,ZhizhangChen,andJiazongZhang, Aencetime-domain methodwithoutthecourantstabilityconditions ,MicrowaveandGuidedWaveLetters, IEEE 9 (1999),no.11,441{443. 107