HIGH-ORDERUNSTAGGEREDCONSTRAINEDTRANSPORTMETHODSFOR MAGNETOHYDRODYNAMICEQUATIONS By QiTang ADISSERTATION Submittedto MichiganStateUniversity inpartialentoftherequirements forthedegreeof AppliedMathematics{DoctorofPhilosophy 2015 ABSTRACT HIGH-ORDERUNSTAGGEREDCONSTRAINEDTRANSPORT METHODSFORMAGNETOHYDRODYNAMICEQUATIONS By QiTang Theidealmagnetohydrodynamic(MHD)equationsareoneofthemostimportantplasma models.Theequationsmodelthedynamicsofaperfectlyconductingquasi-neutralplasma andprovideevolutionequationsforthemacroscopicquantitiesofmass,momentum,and energydensity,aswellasthemagneticMHDhavebeenusedsuccessfullyinmany plasmaphysicsapplicationareas,includinginspaceweatherprediction,astrophysics,as wellasinlaboratoryplasmaapplicationssuchaswsintokamaksandstellarators.Inthis thesis,wefocusonthedevelopmentofhigh-ordernumericalmethodsfortheidealMHD equationsanditsapplications. Inthepartofthethesis,wedevelopaclassofhigh-orderweighted essentiallynon-oscillatory(FD-WENO)schemesforsolvingtheidealMHDequations.In theproposedmethods,wecontroldivergenceerrorsinthemagneticbyusinganovel high-orderconstrainedtransportapproachtosolvethemagneticpotentialequations.The potentialequationsaresolvedusingamoversionoftheFD-WENOschemedeveloped forHamilton{Jacobiequations.Speciallimitersbasedonlresistivityarealsointro- ducedtohelpcontrolunphysicaloscillationsinthemagneticSeveraltwo-dimensional andthree-dimensionalnumericalexamplesarepresentedtodemonstratetheperformance oftheproposedmethod.Numericalresultshaveshownthatwithsuchmethodsweare abletoresolvesolutionstructuresthatareonlyvisibleatmuchhighergridresolutionswith lower-orderschemes. Inthesecondpartofthethesis,wefocusontheproblemsinvolvinglowdensityand lowpressureintheidealMHDsystem.Amaximum-principle-preservinglimiterfor scalarhyperbolicconservationlawsisextendedtoanovelpositivity-preservinglimiterfor theidealMHDequationsinthisportion.TheproposedlimiterisappliedtotheidealMHD schemesproposedinthepart,resultinginahigh-orderpositivity-preservingscheme. Theresultingschemecanachievehigh-orderaccuracy,adiscretedivergence-freecondition andpositivityofthenumericalsolutionsimultaneously.Comparedtotheotherpositivity- preservinglimiterintheliterature,ourlimiterhastheadvantagethatthereisnoextraCFL restrictionfromthelimitingsteps.Numericalexamplesinonedimension,twodimensions andthreedimensionsareprovidedtoverifytheorderofaccuracyonsmoothtestproblems andtoshowtheperformancewhentheproblemsinvolvelowdensityand/orlowpressure. ACKNOWLEDGMENTS Iwouldliketobeginbyexpressingmysinceregratitudetomyadvisor,AndrewChristlieb. Thisthesiswouldnothavebeenpossiblewithouthisexpertise,patience,guidanceand supportoverthepasteyears.Hisinterestinnewscientproblems,hisinsightinthe areaIstudiedandhisinvolvementwithstudentsweregreatlyhelpful.Learningandworking underhisguidancewastrulyanenjoyableexperience. IwouldliketothankYingdaChengforherencouragementandsupportduringmy graduatestudy.Shewasalwaysavailableandprovidedwithgreatideasandkindhelpwhen Iwaslookingforadvices.IamalsogratefultoJamesRossmanithfortheopportunityto workwithhim.Iwanttothankhimforhisintroductiontothetopicsrelatedtothisthesis. Iwouldliketothanktherestofmythesiscommittee:BrianO'Shea,KeithPromislow, andJianliangQian.Iamgratefulforthehelpfuldiscussionswiththemaboutthetopics relatedtothisthesis.Iwouldalsoliketothanksomeother(previous)facultymembersat MSU,includingMattCausley,YamanGu˘clu,WeiGuo,GuanghuiHu,YuanLiu,BenOng, DavidSeal,ZhengfuXu,andXinghuiZhong.Iwouldliketothankmyfellowstudentsat mathdepartmentfortheirfriendshipduringthelasteyears.Thegraduateschoolwould nothavebeensoenjoyableandcolorfulwithoutthem. Lastbutnotleast,Iamgratefulforthesupportofmyfamilyandfriendsduringthe manyyearsasastudent.MydeepestgratitudegoestomywifeZixuan,forbravelyjoining inthisjourneywithme,andforalwayssupportingmeandlovingme.Iamgreatlythankful tomyparents,WeimingTangandMingzhuZhang,fortheirpricelessloveandtrust.Iam alsothankfultomyfamiliesfortheirsupport.Inparticular,Iwouldliketothankmyuncle's family:YanmingZhangandQunShi,aswellasmytwolovelycousins,DannyandSean. iv TheywerealwaystheretohelpwheneverIneededadvicesnomatterinacademiccareeror personallife. v TABLEOFCONTENTS LISTOFTABLES .................................... viii LISTOFFIGURES ................................... ix Chapter1Introduction ............................... 1 1.1Plasma.......................................1 1.2Modelsofplasma.................................2 1.2.1Kineticmodels..............................2 1.2.2Fluidmodels...............................3 1.2.2.1Equationsofmodels...................3 1.2.2.2Magnetohydrodynamicmodel.................5 1.2.2.3Othermodels.......................8 1.3Reviewofpreviouswork.............................9 1.3.1NumericalmethodsforidealMHD...................9 1.3.2Numericalmethodsforhyperbolicconservationlaws..........13 1.3.3Positivity-preservinglimiter.......................14 1.3.4Adaptivemesht........................16 1.4Outlineofthethesis...............................17 Chapter2IdealMHDequations ......................... 19 2.1TheidealMHDequations............................19 2.2Hyperbolicityofthegoverningequations....................20 2.3Magneticpotentialin3D.............................21 2.4Magneticpotentialin2D.............................24 2.5GeneralframeworkofCTalgorithm.......................25 Chapter3FD-WENOwithconstrainedtransportforidealMHD ..... 28 3.1SpatialdiscretizationofidealMHD.......................29 3.2Temporaldiscretization..............................35 3.3Spatialdiscretizationofthemagneticpotentialequation...........36 3.3.1WENOfor1DHamilton{Jacobi.....................37 3.3.22Dmagneticpotentialequation.....................42 3.3.33Dmagneticpotentialequation.....................43 3.4Centraldiscretizationof r A ................49 3.4.1Curlin2D.................................49 3.4.2Curlin3D.................................50 3.4.3Importantproperties...........................51 3.5Numericalresults.................................52 3.5.1SmoothAlfvenwaveproblem......................52 3.5.1.12Dproblem...........................52 vi 3.5.1.23Dproblem...........................53 3.5.22Drotatedshocktubeproblem.....................54 3.5.32DOrszag-Tangvortex..........................57 3.5.4Cloud-shockinteraction.........................60 3.5.4.12Dproblem...........................60 3.5.4.22.5Dproblem..........................62 3.5.4.33Dproblem...........................64 Chapter4Positivity-preservinglimiterforidealMHD ........... 66 4.11Dcase......................................66 4.2Multi-Dcase....................................75 4.3Temporaldiscretization..............................78 4.4Numericalexamples................................80 4.4.1Testcasesin1D..............................80 4.4.1.1Vacuumshocktubetest....................80 4.4.1.2TorsionalAlfvenwavepulse..................81 4.4.2Testcasesinmulti-D...........................82 4.4.2.1SmoothvortextestinMHD..................83 4.4.2.2Smoothvortextestinhydrodynamics............84 4.4.2.3Rotatedvacuumshocktubeproblem.............86 4.4.2.42Dblastproblem........................88 4.4.2.53Dblastproblem........................93 Chapter5Conclusionandfuturework ...................... 96 5.1Conclusion.....................................96 5.2Futurework....................................98 APPENDIX ........................................ 100 BIBLIOGRAPHY .................................... 104 vii LISTOFTABLES Table3.1:Convergencestudyofthe2DAlfvenwaveforCFL=3 : 0.Shown arethe L 2 -errorsand L 1 -errorsattime t =1ofthemagneticld andmagneticpotentialascomputedbytheWENO-CT2Dschemeat variousgridresolutions..........................53 Table3.2:Convergencestudyofthe3DAlfvenwaveforCFL=3 : 0.Shownare the L 2 -errorsand L 1 -errorsattime t =1ofallthemagneticand magneticpotentialcomponentsascomputedbytheWENO-CT3D schemeonvariousgridresolutions....................54 Table4.1:Accuracytestofthe2DvortexevolutioninMHD.Shownarethe L 1 -errorsand L 1 -errorsattime t =0 : 05ofthedensityascomputed bythepositivity-preservingWENO-CT2Dschemeatvariousgrid resolutions................................85 Table4.2:Accuracytestofthe2Dvortexevolutioninhydrodynamics.Shown arethe L 1 -errorsand L 1 -errorsattime t =0 : 05ofthedensityas computedbythepositivity-preservingWENO-CT2Dschemewith CTstepsturnedThesolutionsconvergeataccuracy.86 Table4.3:Comparisonsoftschemessolvingthe2Dblastproblem.The columnof P indicateswhetherthenumericalsolutionsremainposi- tiveinthesimulations.Thecolumnof S indicateswhetherthesim- ulationsrunstablyto t =0 : 01.Inordertomakeafaircomparison, thepositivityofdensityandpressureisonlycheckedateachtime step t n ..................................91 viii LISTOFFIGURES Figure3.1:1Dadvectionequation.Showninthesepanelsare(a)thesolution obtainedbytheWENO-HCLschemewith(b)itsderivative,and(c) thesolutionobtainedbytheWENO-HJschemewith(d)itsderiva- tive. Forinterpretationofthereferencestocolorinthisandallother es,thereaderisreferredtotheelectronicversionofthisdisser- tation. ..................................41 Figure3.2:Therotatedshocktubeproblem.Showninthesepanelsare ˆ of (a)thebaseWENOschemeand(b)theWENO-CT2Dscheme.30 equallyspacedcontoursareshownforeachgraphintheranges ˆ 2 [0 : 1795 ; 1].................................55 Figure3.3:Therotatedshocktubeproblem.Showninthesepanelsare(a)a1D cutalong y =0of ˆ ascomputedwiththebaseWENOscheme,(b) azoomedinviewofthesameplot,(c)a1Dcutalong y =0of ˆ as computedwiththeWENO-CT2Dscheme,and(d)azoomedinview ofthesameplot.Thesolidlineisahighlyresolved1Dsolution...56 Figure3.4:Therotatedshocktubeproblem.Showninthesepanelsare1Dcuts along y =0ofthemagnetic(a)perpendicular( B ? )and(b)par- allel( B k )totheshockinterfaceascomputedwiththebaseWENO schemeand1Dcutsalong y =0ofthemagnetic(c)perpendic- ular( B ? )and(d)parallel( B k )totheshockinterfaceascomputed withtheWENO-CT2Ddscheme.Thesolidlineisahighlyresolved 1Dsolution................................57 Figure3.5:TheOrszag-Tangvortexproblem.Showninthesepanelsare ˆ at times(a) t =0 : 5,(b) t =2,(c) t =3,and(d) t =4ascomputed withtheWENO-CT2Dschemeona192 192mesh.15equally spacedcontourswereusedforeachplot................58 Figure3.6:TheOrszag-Tangvortexproblem.Showninthesepanelsare(a) thethermalpressureascomputedwiththeWENO-CT2Dschemeat t =3ona192 192meshand(b)asliceofthethermalpressureat y =0 : 625 ˇ ................................59 Figure3.7:The2Dcloud-shockinteractionproblem.Showninthesepanelsare schlierenplotsattime t =0 : 06of(a)ln( ˆ )and(b)thenormof B .ThesolutionwasobtainedusingtheWENO-CT2Dschemeona 256 256mesh.............................61 ix Figure3.8:The2Dcloud-shockinteractionproblem.Showninthesepanelsare B z attime t =0 : 06ona256 256meshassolvedwith(a)the WENO-CT2Dschemeand(b)the2.5DversionoftheWENO-CT3D scheme.25equallyspacedcontourswereusedineachofthesepanels.63 Figure3.9:The3Dcloud-shockinteractionproblem.Showninthisplotisln( ˆ ) ascomputedwiththeWENO-CT3Dschemeona128 128 128 mesh...................................65 Figure4.1:Thevacuumshocktubeproblem.Showninthesepanelsareplotsat time t =0 : 1of(a) ˆ and(b)thethermalpressure.Thebluecircleis asolutionsolvedonameshwith N =200.Thesolidlineisahighly resolvedsolution.............................81 Figure4.2:ThetorsionalAlfvenwavepulse.Showninthesepanelsareplotsat time t =0 : 156of(a)theenergyand(b)thethermalpressure.The solutionwasobtainedonameshwith N =800............82 Figure4.3:ThetorsionalAlfvenwavepulse.Showninthesepanelsareplotsat time t =0 : 156of(a) u y ,(b) u z ,(c) B y and(d) B z .Thesolution wasobtainedonameshwith N =800.................83 Figure4.4:Rotatedvacuumshocktubeproblem.Showninthesepanelsareplots attime t =0 : 1of(a) ˆ ,(b) ˆ cutaty=0,(c)thethermalpressure and(d)thethermalpressurecutaty=0.40equallyspacedcontours areusedforthecontourplots.Thesolidlinesin(b)and(d)are1D highlyresolvedsolutions.Thesolutionwasobtainedona240 100 mesh...................................87 Figure4.5:2Dblastproblem.Showninthesepanelsareplotsattime t =0 : 01of (a) ˆ ,(b)thethermalpressure,(c)thenormof u and(d)themagnetic pressure.40equallyspacedcontoursareusedforeachplot.The solutionwasobtainedona256 256meshbypositivity-preserving WENO-CTscheme...........................89 Figure4.6:2Dblastproblem.Showninthesepanelsareplotsattime t =0 : 01 of(a) ˆ ,(b)thethermalpressure,(c)thenormof u and(d)the magneticpressure.40equallyspacedcontoursareusedforeachplot. Thesolutionwasobtainedona256 256meshbytheWENO-HCL scheme..................................90 x Figure4.7:ComparisonsbetweenWENO-HCLschemeandpositivity-preserving WENO-CTschemeforthe2Dblastproblem.Showninthesepanels areplotsof(a)thedivergenceerrorand(b)therelativeerrorof thetotalenergyinthewholedomain.The x -axisdenotesthetime t 2 [0 ; 0 : 01]................................92 Figure4.8:3Dblastproblem.Showninthesepanelsare3Dplotsattime t =0 : 01 of(a) ˆ and(b)thethermalpressure.Thesolutionwasobtainedon a150 150 150mesh.........................93 Figure4.9:3Dblastproblem.Showninthesepanelsareplotsattime t =0 : 01 andcutat z =0of(a) ˆ ,(b)thethermalpressure,(c)thenormof u and(d)themagneticpressure.40equallyspacedcontoursareused foreachplot.Thesolutionwasobtainedona150 150 150mesh.94 xi Chapter1 Introduction 1.1Plasma Inphysics,matterintheuniverseisusuallyintofourstates: solid , liquid , gas and plasma .Thefundamentaldistinctionsamongthosefourstatesarebasedonof thestrengthofthebondsbetweentheparticlesofthestates.Itisalsoequivalenttosay themolecularinterrelationshipsmatter.Astateiscalled solid ifintermolecularattractions keepthemoleculesinspatialrelationshipssuchthatthebindingforcesareverystrong. Liquid isthestateinwhichtheweakattractionskeepmoleculesclose,butnotin relationships.Whenmoleculesarerelativelyseparatedandintermolecularattractionsare essentiallyabsent,moleculesareinthegaseousstate.However, plasma isanionizedgasthat occursatveryhightemperatures.Theatomsormoleculesinthisheatedgashavesomuch kineticenergythattheyareabletoovercomethebindingenergyoftheelectrons.Although thetransitionfromagastoaplasmaisnotaphasetransitioninthethermodynamicsense, thedistinctintermolecularforcesbytheheatingprocessleadstothedistinctpropertiesof plasma.Asaresult,plasmaistreatedasthefourthstateofmatter.OnEarth,thecommon statesofmatteraresolid,liquidandgas,butmostofthematteroftheuniverseisfound intheplasmaform.Forinstance,starsaremostlymadeofplasma.Inthisthesis,weare generallyinterestedinthemodelsofplasma. 1 1.2Modelsofplasma Therearemainlytwodistincttypesofmodelsforplasma: kineticmodels and models . Bothofthemdescribetheparticlelocationsandvelocitiesandtheelectromagnetic intheplasma.Thekineticmodelsdescribeplasmabyadistributionfunctionofparticles ateachpointinthe phasespace ,whilethemodelsusemacroscopicquantitiessuchas massdensityandpressuretodescribeplasma.Ingeneral,kineticmodelsareexpectedtobe moreaccuratethanmodels.Howeveritisalsomorecomputationallyexpensivedueto itshigherdimension.Inthissectionwereviewseveralkineticandmodels. 1.2.1Kineticmodels Ingeneral,kineticmodelsrefertothemethodsevolvingsomerepresentationofparticle positionsandvelocities.Oneofthemostfundamentalequationstodescribeplasmais the Boltzmannequation ,whichdescribestheevolutionofaprobabilitydensityfunction f ( t; x ; v ).Here f ( t; x ; v )isadistributionfunctionofplasmaspecies inthephasespace ( x ; v )atagiventime t .Mathematically,itiswrittenasfollows: @f @t + v r f + a r v f = @f @t coll (1.1) Inthecasewherelongrangeinteractionsbetweenchargedparticlesaredescribedwith Maxwell'sequations,theacceleration a istypicallydependentonlyonthe Lorentzforce by theelectromagnetici.e., a = F m = q m ( E + v B ) ; (1.2) 2 where m representsthemassofspecies and q representsthechargeof .Theterm @f @t coll in(1.1)representsageneralcollisionoperator,whichdescribesthechangeof plasmaduetoparticlecollisions.Inparticular,the Vlasovequation isthecasewhenthe collisionsbetweenparticlesareneglected,whichcanberesultedundercertainscalingofthe theBoltzmannequation.Inotherwords,weassume @f @t coll =0resultingintheVlasov equation. 1.2.2Fluidmodels Aswementioned,thekineticmodelsarenumericallyexpensivetosolveduetothehigh dimensionofthephasespace.Asaresult,peoplesimplifythehigh-dimensionalBoltzmann equationbyevaluatingthemomentsinthevelocityspace( v -space).Withaproperclosure forhighermomentsusedintheprocessofmomentevaluation,theresultingsystemsbecome themodelswhoseunknownsonlydependontheindependentvariablesofthetime t andthespace x .Inthissubsectionwepresentageneralproceduretoobtain modelsandandalsoreviewseveralimportantmodelsofmodels. 1.2.2.1Equationsofmodels Ingeneral,theunknownsinmodelsarebyintegralofmomentsof v overthe velocityspace.Theircorrespondingequationscanbeobtainedbyintegralofproductof theBoltzmannequation(1.1)withthose v moments.Forinstance,multiplyingtheabove equationbysomefunction ˜ ( v )andintegratingoverthewhole v -space,weget, Z ˜ @f @t d v + Z ˜ v r f d v + Z ˜ a r v f d v = Z ˜ @f @t coll d v : (1.3) 3 Rearrangingtheaboveequationas: @ @t Z ˜f d v + r Z ˜ v f d v + Z ˜ ( a r v ) f d v = @ @t Z ˜ ( f ) coll d v (1.4) Ifwedethenotation <˜> astheaverageover v -space <˜> = Z ˜f d v ; then(1.4)willbecomeageneralizedtransportequation: @ @t <˜> + r ( <˜ v > ) < ( a r v ) ˜> = @ @t ( <˜> ) coll (1.5) Asanexample,ifwelet ˜ bethemomentsof m , m v and 1 2 m k v k 2 respectively,weobtain threetransportequationsfromequation(1.5): @ˆ @t + r ( ˆ u )= @ˆ @t coll (1.6) @ˆ u @t + r ( ˆ u u )+ r ( P ) < F > = A (1.7) @ @t + r ( u )+ r ( P u ) n q u E = M (1.8) Heretheunknownsarethemassdensity ˆ = n m = m < 1 > ,themomentum density ˆ u = ˆ < v > andtheenergydensity bythepressure p andthe kineticenergy, = p 1 + 1 2 ˆ k u k 2 : 4 with beingtheidealgasconstant.Othertermsintheequationsarethepressuretensor P = m < ( v u ) ( v u ) > andthecollisionterms, A = m Z v @f @t coll d v ;M = 1 2 m Z k v k 2 @f @t coll d v : Ingeneral,thisprocesscanbeextendedtoarbitrarynumberofmomentsof v .Toobtain aclosedsystem,wehavetocutthesystemforacertainnumberofmoments.However,the transportequationofthehighestmomentintheresultingsystemwillrequireanevenhigher momentinitsfunction.Thisleadstotheso-called closureproblem ,inwhichtheof thehighestmomentneedstobeapproximatedbyallthemomentsintheresultingsystem. Thisapproximationhastobereasonableinthesensethatthesystemissolvableandthe solutionhasphysicalmeaning. 1.2.2.2Magnetohydrodynamicmodel Plasmagenerallyinvolvesmultiplespeciesofparticles.Onespecieswillleadtoonesystem followingtheprocedurediscussedinSection1.2.2.1.Ingeneral,thosesystemseach otherthroughthecollisionoperatorandtheelectromagneticThemagnetohydrody- namic(MHD)equationsfurtherassumethereisonlyoneinthesystemduetothescale oftheconsideredproblem.Inotherwords,theunknownsintheMHDequationsarethe summationwithrespecttooftheunknownsfromSection1.2.2.1,forinstance,thetwo 5 unknownsare massdensity: ˆ ( t; x )= X n m ; momentumdensity: ˆ u ( t; x )= X n m u : Withtheadditionalassumptiontosimplifythepressuretensor P (thisassumptionisac- tuallytheclosurecondition),andnocollisionsontherighthandside,wecanobtainthree conservationequationsforthemass,momentumandkineticenergydensitiesbysimplysum- mingthecorrespondingequationforeachspecies.Thosethreeequationsareexactlythe threeconservationlawsintheidealMHDequations. ThelastequationoftheMHDequationscomesfromFaraday'slaw, @ B @t + r E =0 : tfromthefullMaxwellequationcase,theelectricintheMHDequationis directlydeterminedbytheOhm'slaw. Forinstance,wecouldconsiderthecasewhenthereareonlytwospecies,ionsandelec- trons,intheplasma.FollowingtheprocessinSection1.2.2.1,wecouldobtaintwosystems oftransportequationsforbothspecies.Wefurtherassumetheunknownssuchasthecurrent density J andthechargedensity ˆ c , J = e ( n i u i n e u e ) ; ˆ c = e ( n i n e ) : 6 Additionally,thereisanotherimportantparameter.The Debyelength D isdbythe typicalnumberdensity n 0 andthetemperaturescale T 0 , D = 0 k B T 0 e 2 n 0 1 2 : IntheMHDlimit,thespatialscaleoftheproblemismuchlargerthan D .FromGauss's law,itisequivalenttoassumetheplasmais quasi-neutral ,i.e., ˆ c =0or n i = n e = n . Bytakingtheoftwomomentumtransportequations,wecouldobtainthe generalizedOhm'slaw asfollows, @ J @t + r e m i p i e m e p e I + uJ + Ju JJ en m i m e m i + m e ˆe 2 m i m e ( E + u B )+ e ( m i m e ) m i m e J B =0 ; (1.9) whichcanberewrittenas, E + u B = m i m e ˆe 2 @ J @t + r e m i p i e m e p e I + uJ + Ju JJ en m i m e m i + m e + m i m e ˆe J B : (1.10) IntheidealMHDlimit,theratio m e m i isassumedtobeverysmall.Ifwetakethelimit m e m i =0,theequation(1.10)canbeas, E + u B = m e ne 2 @ J @t 1 ne r p e + m e ne 2 r uJ + Ju JJ en + 1 ne J B : (1.11) IntheidealMHDlimit,therighthandsideoftheequation(1.11)isignoredduetothe spatialscaleoftheproblem.Onewaytojustifythatistoassumethespatialscalelengthis 7 muchgreaterthantheions'inertiallength( c ! i = r m i 0 ne 2 ).Thisresultsinthe idealOhm's law , E + u B =0 : (1.12) Ifothertermsin(1.9)areinvolvedintheproblem,thesystemisgenerallycalled extended MHD .Forexample,if J B (Hallterm)isincluded,thewholesystemiscalled HallMHD .In thisthesis,weonlyfocusontheidealMHDequationsandtheextensiontootherextended MHDequationswillbepartofourfuturework. 1.2.2.3Othermodels TheMHDequationsarestilllimitedintheirabilitytofullydescribeaplasma.Forinstance, intheidealMHDequations,themagneticlinesmovewiththeandthuswillnot changeitstopologyastimeevolves.However,duetothemomentumoftheelectrons,the themagneticlineinmanyproblemscanbreakandreconnect.Sincetheelectronsare masslessintheidealMHD,thisectisapparentlybeyondthelimitoftheidealMHD equations.Sotostudytheproblemsofthisnature,suchasthesemagneticreconnection problems,otheruidmodelshavebeenproposedandstudied. AnaturalgeneralizationoftheMHDequationsisthe model ,wheretheelectrons andionsaremodeledbytwosetsofvariables.Dependingontheclosureoftheequations, thetwmodelshaveseveralvariations.Ifthesystemonlyevolvestheenergyorentropy equationsseparatelywhilestillusingmassandmomentumdensityequationsand Ohm'slaw,thewholesystemisoftencalledthe MHDequations .Howeveramore generalcaseistousethefullMaxwellequation,inwhichcasetheelectronsandionsare modeledbycompletelyseparatevariables.Fortwmodels,Johnsonhasagood 8 summaryofentmodelsinhisthesis[41]. AnotherapproachistogeneralizetheprocedureoftheBoltzmannequationto models.Insteadofevolvingonlythreeconservationvariables,themodelusesmorehigher momentof v .Thisisoftencalled momentmethods .Themotivationofthemomentmethods isthatwithmoremomentsintroduced,theresultingsystemwillbecomeabetterapproxima- tiontotheBoltzmannequationthantheMHDsystem.Itisadvantageousthatthesystem stillhasthesamedimensionastheMHDequations.Grad's13-momentsystem[32]isone oftheearlyworksinthisarea.Howevertherearemanyinthesemethods.For instance,itisgenerallyhardtoobtainahyperbolicsystembyusingtheclosuretoclosethe systemforhighmomentsystems.Recently,therehavebeensomeincreasinginterestinthis approachbecausepeopleproposednovelapproachestoclosethesystem,e.g.,[46,72,74]. Butthisapproachstillhasalongwaytogobeforeitreachesalevelofpracticalapplication. 1.3Reviewofpreviouswork Inthissection,wereviewthepreviousworkrelatedtonumericalmethodsfortheMHD equations.Wemainlyfocusonfourtopics:(1)numericalmethodsforidealMHD,(2) adaptivemesht(AMR),(3)numericalmethodsforthehyperbolicconservation laws,and(4)positivity-preservinglimiters. 1.3.1NumericalmethodsforidealMHD Mathematically,theidealMHDequationsareasetofnonlinearhyperbolicconservationlaws withtheadditionalrestrictionthatthemagneticmustremaindivergence-freeforall time.Infact,atthecontinuumlevel,iftheinitialmagneticisdivergence-free,then 9 theidealMHDequationspropagatethisconditionforwardforalltime.Unfortunately,most standardnumericaldiscretizationsbasedonshock-capturingmethods(e.g.,volume, weightedessentiallynon-oscillatory,discontinuousGalerkin)donotpropagateadiscrete versionofthedivergence-freeconditionforwardintime;andfurthermore,therearesome observationsintheliteraturethatthefailuretoguarantee r B =0toerror canleadtosomeunphysicalsolutions.Forinstance,BrackbillandBarnes[16]the unphysicalwduetotheerrorfrom r B eventuallyeventerminatesthecomputation. AnotherexampleistheRotorproblemin[10]thatshowsthesolutionswillbespurious withoutanydivergencecleaningstep. ThemainchallengeinnumericallysimulatingtheidealMHDsystemisthereforeto augmentexistingschemessothattheysatisfyadivergence-freeconditionofthemagnetic insomeway.Roughlyspeaking,therearefourkindsofapproachesthathavebeenproposed intheliterature:(1) the8-waveformulation [56,57],(2) projectionmethods [6,16,75,83], (3) hyperbolicdivergence-cleaningmethods [26],and(4) constrainedtransport(CT)methods [3,10,22,25,28,29,34,51,52,60,62,71,75,73]. Themethodin[16]isprobablytherstapproachdealingwiththisconstraintusing theclassicalprojectionmethod.Theyused Hodgedecomposition toprojecttheinaccurate magnetictothedivergence-freesubspacebysolvingaPoissonequation.However, thisprojectionmethodiscomputationallyexpensiveandgenerallytoextendtothe AMRframeworks.Therearetwofactorsdrivinghere,theneedarisestosolveaPoisson equationeachtimestep,second,solvingthePoissonequationonanadaptivemeshcan betricky.Thistreatmentisstillverypopularbecauseofitsaccuracyandxibility.For instance,Zachary etal. [83]splitstheidealMHDEquationsintoaconservativepartand anotheradvectivepart,treatedasasource.AnewRiemannsolverwasintroducedtoget 10 ridofmostofunphysicalstatessuchasnegativepressuresanddensities,andtheprojection methodwasusedtomakethemagneticsatisfytheconstraint.Thesecondapproach inventedbyPowellin[56]istotreattheidealMHDEquationstlyt,by addinganadditionalsourceterm(itisproportionaltothedivergenceofthemagnetic totheMHDequations.Consequently,theoriginalidealMHDequationswiththeconstraint ischangedintoa8 8hyperbolicsystemwithasourcetermandnoconstraint.Therefore, thissystemiseasiertosolveandtheresulting 8-wavescheme isveryrobustandextension toanAMRframeworkismucheasierforthismethod.Thedrawbackofthisapproachis alsoveryclear.Thedivergenceisnotenforcedtomachineprecisionandthenewsystem isnotconservativeanymoreduetotheadditionalsourceterm.Ithasbeenobservedin [75]thatthisnon-conservationdrawbackcanproduceincorrectjumpsinarotatedshock tubeproblem.Anotherclassofmethodsisthe hyperbolicdivergence-cleaningmethod proposedbyDedner etal. [26],whichhasasimilarideaastheprojectionmethod.Instead ofsolvinganellipticequationintheprojectionmethod,amixedhyperbolicandparabolic equationisintroducedandsolvedforthedivergenceerrorofmagneticThemethod dampsthedivergenceerrorsaway,insteadofgettinganexactdivergence-freemagnetic Theadvantageofthismethodisthatitisfullyexplicit,thus,tandfast,butthe disadvantageisthatithastwotunableparameters:thespeedofpropagationoftheerror andtherateatwhichthedivergenceerrorisdamped. InthisthesiswefocusontheCTmethodologyforproducingamagneticthat adiscretedivergence-freecondition.TheCTmethodwasoriginallyintroducedby EvansandHawley[28],and,intheirformulation,staggeredelectricandmagneticare usedtocreateappropriatemimetic(FD)operatorsthatultimatelyleadto anexactlydivergence-freemagneticTheirconstrainedtransportframeworkcanbe 11 thoughtasamoofthepopularYeescheme[82]fromelectromagnetics. SincetheintroductionoftheCTmethodology,therehavebeenmanymoand extensions,especiallyinthecontextofhigh-resolutionshock-capturingschemes.DeVore[27] developedatransportimplementationoftheconstrainedtransportapproach. BalsaraandSpicer[10],DaiandWoodward[25],andRyu etal. [62],alldevelopedvarious strategiesforconstructingtheelectricviaOhm'slawintheCTframework.Londrillo andZanna[51,52]proposedahigh-orderversionoftheEvansandHawleyapproach.De Sterck[71]developedasimilarCTmethodonunstructuredtriangulargrids.Balsara[1] developedfromtheCTframeworkanAMRschemethatalsoincludedagloballydivergence- freemagneticldreconstruction.Thereisacarefuldescriptionandcomparisonofseveralof thesemethodsinthearticleofoth[75],inwhichhealsoshowedthatastaggeredmagnetic isnotnecessary,andthenintroducedseveralunstaggeredCTmethods. Inrecentyears,unstaggeredCTmethodshaveattractedconsiderableinterestduetotheir easeofimplementationandapplicabilitytoAMRstrategies.Forinstance,FeyandTorrilhon [29]presentedawaytopreservedivergence-freeconditionthroughanunstaggeredupwind scheme.Rossmanith[60]developedanunstaggeredCTmethodforthe2DMHDequations onCartesiangridsbasedonthewave-propagationmethod[45].Helzel etal. [34,35]extended thisunstaggeredCTmethodtothe3DMHDequationsandtomappedgrids. Inadditiontotheabovementionedpapers,severalotherhigh-ordermethodshavebeen proposedinrecentyearsfortheidealMHDequationsusingavarietyofdiscretizationtech- niques.Balsara[4]developedaweightedessentiallynon-oscillatory(WENO)methodfor idealMHDusingastaggeredmagnetictoreconstructagloballydivergence-freemag- neticBalsara etal. [7,8]developedaclassofhigh-orderADER-WENOschemes,again usingastaggeredmagnetictoreconstructagloballydivergence-freemagneticLi 12 etal. [49]andCheng etal. [19]introducedaclassofcentraldiscontinuousGalerkinschemes thatevolvestheMHDequationsonaprimalaswellasadualmesh,andbyintertwiningthese twoupdates,theyshowedthatagloballydivergence-freemagneticcouldbeobtained. Kawai[42]devlopedahigh-ordermethodwithresistivity,where theeoperatorswerespconstructedtoguaranteethatanappropriate ofthemagneticdivergenceispropagatedforwardintimebythenumerical scheme. 1.3.2Numericalmethodsforhyperbolicconservationlaws Inthisthesis,weusethe methodoflines approachtosolvethehyperbolicconservation laws.Theideaofthisapproachistodiscretizethetime-involvedPDEsonlyinspace togenerateasemi-discretizedsystemofODEs,andthentoapplyanODEsolvertothe resultingsystem.Manynumericalmethodshavebeendevelopedforspatialdiscretization oftheconservationlawovertherecentdecades,suchasthediscontinuousGalerkinmethod [23],thevrenceENOschemes[33,69,70],andv WENOschemes[39,50].Amongvariousmethods,WENOschemesareshown tobeveryrobustandtespeciallywhensolutionsmaycontaindiscontinuities,sharp gradientregionsandothercomplicatedsolutionstructures.Inthisthesis,weuseWENO methodtoserviceasafundamentalsolverfortheidealMHDequations. Forthetimeintegrators,thenon-oscillatorypropertyisdesiredfortheproblemsofthe conservationlaws,suchasEulerequationsandMHDequations.Thus,strongstability preservingRunge{Kutta(SSP-RK)methods(alsocalledtotal-variation-diminishinginthe earlierliterature)areverypopularforthoseproblems.SSPmethodsarehigher-ordermeth- odsthatpreservethestrongstabilitypropertiesofEulertimesteppingforthe 13 spatialdiscretization.TheoptimalthirdorderSSP-RKmethod(SSP-RK3)isgivenby[69] q (1) = q n + L ( q n ) q (2) = 3 4 q n + 1 4 q (1) + 1 4 L ( q (1) ) q n +1 = 1 3 q n + 2 3 q (2) + 2 3 L ( q (2) ) (1.13) InChapter4,weusetheaboveSSP-RK3methodasourtimeintegratortodesignapositivity- preservinglimiter.InChapter3,weusealow-storagefourth-orderSSP-RKmethodtosolve theproblems,whichispresentedinthework[43]. 1.3.3Positivity-preservinglimiter Anothermajorfocusofthisthesisisthedesignofhigh-orderschemesthatpreservethepos- itivityofthedensityandpressureoftheMHDsystem.Evenwithdivergence-freemethods, negativedensityor/andpressurecanstillbeobservedinnumericalsimulations,suchasthose simulatingthelow- plasma.Thisnegativequantitycanleadtoacomplexwavespeedthat breaksthehyperbolicityofthesystemandcausesthenumericalsimulationstobreakdown. Alotofhasbeendedicatedtoaddressingthisissueintheliterature.Forinstance,Bal- saraandSpicer[9]proposedastrategytomaintainthepositivityofpressurebyswitchingthe Riemannsolversbasedontwavesituations.Janhunen[37]designedanewRiemann solverforthemoidealMHDequationsanddemonstrateditspositivity-preserving propertynumerically.In[76],aconservativesecond-orderMUSCL-Hancockschemewas showntobepositivity-preservingforthe1DidealMHDequationsandtheextensionto multidimensional(multi-D)caseswasconstructedbasedonsimilarideasasPowell's8-wave formulation[56,57].Balsara[5]developedahigh-orderpositivity-preservingschemeforideal 14 MHDthroughlimitinghigh-ordernumericalsolutionsbyaconservativeboundedsolution. AnotherclassofimportantmethodsfortheidealMHDequationsisdiscontinuousGalerkin (DG)methods[47,48,49,61,81].Recently,Cheng etal. proposedpositivity-preserving DGandcentralDGmethodsfortheidealMHDequations[19],inwhichtheygeneralized ZhangandShu'spositivity-preservinglimitersforthecompressibleEulerequations[84].In [19],itwasprovedthattheLax-Friedrichsschemeispositivity-preservingforthe 1DMHDundertherestrictionCFL 0 : 5.Thisschemealsoservesasthebuilding blockforthepositivity-preservinglimiterinthisthesis. BesidestheaforementionedworkforMHDequations,severalhigh-orderpositivity-preserving schemeshavebeendevelopedrecentlyforcompressibleEulerequations.ZhangandShu developedarbitrary-orderpositivity-preservingvolumeWENOandDGmethodsby limitingtheunderlyingpolynomialsaroundcellaverages[84].Alimiterwas proposedbyHu etal. [36]forWENOschemestomaintainpositivityof densityandpressureforthecompressibleEulerequations.Recently,Xuhasproposeda maximum-principle-preservinglimiterforWENOschemesforthescalar conservationlawsin[79].Thislimiterhaslaterbeenextendedtoapositivity-preserving limiterforteWENOschemes[78].Intherecentworkof[20],wegeneralized thislimitertohighordervolumeWENOschemesontriangularmeshesforthescalar conservationlawsandEulerequations.Inanotherrecentwork[65],wesuccessfullyapplied theabovepositivity-preservinglimitertothesingle-stepmethodforEulerequations[63]. Theapproachdevelopedin[20,65,78]isnovelbecausetheparametrizedlimitermaintains theaccuracyofthebaseschemewithoutsacritheCFLexcessivelyandismoret thantheotherversionsofthelimiters[36,84]. 15 1.3.4Adaptivemesht TheAMRalgorithmwasdevelopedforhyperbolicconservationlawsin[15].Theoriginal ideaistothemesharoundcomplicatedsolutionsstructure,suchasshocks,resultingin abetterresolutionanddecreaseofcomputationalIn[15],afundamentalframework ofpatch-basedAMRhasbeenestablished.Theoriginalmethodstartsfromauniformrect- angularmesh,usingRichardsonextrapolationtoestimatelocalspacialerrorstodetermine thementregions.Aftertheregionsarethealgorithmsolvestheoriginalhyper- bolicequationslocallywithappropriatelocalboundaryconditions.Sincethealgorithmalso timestepslocally,thecomputationisverytbecausetheCFLofthecoarse meshesisnotrestrictedbytheCFLofthemeshes. Sincethen,theAMRalgorithmhasbeenverypopularinalargevarietyofareassuch asspacephysics,atmosphericmodelingandaerospaceengineering.Accompanywiththe applications,therearealsomanyondevelopingthetheoriesofAMRalgorithms.In [13],BergerandColellainvestigatedtheglobalconservationofthesolutioninAMR,achieved bythemoonoftheupdatingontheinterfacebetweenandcoarsegrids.Later, Bell etal. [12]extendedthisideato3Dhyperbolicconservationlaws.In[14],thispatched- basedAMRwasalsosuccessfullycombinedwiththevolumewave-propagationmethod thatwasproposedbyLeVequein[45]. TherehavebeenalotofstudiesofAMRalgorithminEulerEquationsandMHDEqua- tions.ApackagecalledGerris[55]isdesignedforincompressibleEulerEquation.The algorithmusesthefullythreadedtreedatastructureof[44]andinventedamulti-gridPois- sonsolvertoimplementtheprojectionmethodbasedonHodgedecomposition.ForMHD equationssimulations,Rossmanith[59]wasabletoextendtheunstaggeredCTmethodto 16 AMRCLAWframeworkfor2DidealMHDequations.AnotherextensionoftheCTmeth- odstoAMRframeworkwasproposedbyFromang etal. [30],inwhichtheyusetheAMR framework,MUSCLschemeandthestandardCTmethodtoevolvetheinductionequations intime.In[2],theCTmethodfrom[10]wasalsoextendedtoAMRhierarchybyBalsara. 1.4Outlineofthethesis Inthework,wemainlyconsidertwonumericalultieswhensolvingtheidealMHD equations.First,wewouldliketosolvetheidealMHDequationinsuchawayastomaintain thedivergence-freepropertyofthemagneticSecond,wewouldliketomaintainthe positivityofbothdensityandpressurewhenthemethodisusedtomodelanyplasma problem. InChapter2weintroducethemathematicalformoftheidealMHDequations.We alsopresentitshyperbolicpropertybydiscussingitseigenvalues.Wefurtherintroducethe generalconstrainedtransportmethodandthe2Dand3Dversionsofthemagneticpotential equationsusedinthiswork. InChapter3weproposeaclassofnovelerenceschemesforthe2Dand3Dmag- neticpotentialequations.TheproposedschemesarecoupledwithregularWENOschemefor hyperbolicconservationlawstoobtainafourth-ordernumericalschemesfor theidealMHDequations.Theschemessuccessfullycontroltheoscillationsinthesolutions thankstotheWENOapproachandnovelschemesweproposedforthemagneticpotential equations.Theresultingschemesaretestedwithseveraltestproblems,including2Dand 3DsmoothAlfvenwaveproblems,2Drotatedshocktubeproblem,Orszag-Tangvortex,and 2D,2.5Dand3Dcloud-shockinteractionproblems. 17 InChapter4weproposeapositivity-preservinglimiterfortheschemesweproposedin Chapter3.ThelimiterisanaturalextensionofasimilarlimiterfortheEulerequations. Thelimitergenerallyhastwosteps:itmaintainspositivityofthedensitybyalinear process,andthenitmodifythetokeepthepressurepositivebyanonlinearprocess. AftertheproposedlimiterisappliedtotheconstrainedtransportschemesinChapter3,the resultingschemecansolvetheproblemswithlowdensityandpressurewithaCFLnumber of0 : 5.Wefurthertestitwithseveralproblems,including1Dand2Dvacuumshocktube problems,torsionalAlfvenwavepulseproblem,smoothvortexproblems,and2Dand3D blastproblems. 18 Chapter2 IdealMHDequations 2.1TheidealMHDequations TheidealMHDequationsinconservationformcanbewrittenas @ @t 2 6 6 6 6 6 6 6 6 6 4 ˆ ˆ u E B 3 7 7 7 7 7 7 7 7 7 5 + r 2 6 6 6 6 6 6 6 6 6 4 ˆ u ˆ u u +( p + 1 2 k B k 2 ) I B B u ( E + p + 1 2 k B k 2 ) B ( u B ) u B B u 3 7 7 7 7 7 7 7 7 7 5 =0 ; (2.1) r B =0 ; (2.2) where ˆ , ˆ u ,and E arethetotalmass,momentumandenergydensitiesofthesystem, B is themagneticand p isthehydrodynamicpressure.Thetotalenergydensityisgivenby E = p 1 + 1 2 ˆ k u k 2 + 1 2 k B k 2 ; (2.3) where =5 = 3istheidealgasconstant.Here kk isusedtodenotetheEuclideanvector norm.AcompletederivationoftheMHDsystem(2.1){(2.2)canbefoundinmanystandard plasmaphysicstextbooks(e.g.,pages165{190of[54]). Thischapterpreviouslyappearedin[22]:A.J.Christlieb,J.A.RossmanithandQ.Tang.Finitedierence weightedessentiallynon-oscillatoryschemeswithconstrainedtransportforidealmagnetohydrodynamics. J. Comput.Phys., 268:302{325,2014. 19 2.2Hyperbolicityofthegoverningequations Equations(2.1),alongwiththeequationofstate(2.3),formasystemofhyperbolicconser- vationlaws: @q @t + r F ( q )=0 ; (2.4) where q =( ˆ;ˆ u ; E ; B )aretheconservedvariablesand F isthetensor(see(2.1)). Undertheassumptionofpositivepressure( p> 0)anddensity( ˆ> 0),theJacobian insomearbitrarydirection n ( k n k =1), A ( q ; n ):= n @ F @q ,isadiagonalizablematrixwith realeigenvalues.Inparticular,theeigenvaluesoftheJacobianmatrixinsomearbitrary direction n ( k n k =1)canbewrittenasfollows: 1 ; 8 = u n c f :fastmagnetosonicwaves,(2.5) 2 ; 7 = u n c a :Alfvenwaves,(2.6) 3 ; 6 = u n c s :slowmagnetosonicwaves,(2.7) 4 = u n :entropywave,(2.8) 5 = u n :divergencewave,(2.9) 20 where a r p ˆ ; (2.10) c a s ( B n ) 2 ˆ ; (2.11) c f 8 < : 1 2 2 4 a 2 + k B k 2 ˆ + s a 2 + k B k 2 ˆ 2 4 a 2 ( B n ) 2 ˆ 3 5 9 = ; 1 2 ; (2.12) c s 8 < : 1 2 2 4 a 2 + k B k 2 ˆ s a 2 + k B k 2 ˆ 2 4 a 2 ( B n ) 2 ˆ 3 5 9 = ; 1 2 : (2.13) Theeighteigenvaluesarewell-orderedinthesensethat 1 2 3 4 5 6 7 8 : (2.14) 2.3Magneticpotentialin3D Althoughtherearemanyavailablenumericalmethodsforsolvinghyperbolicsystems(e.g., volume,WENO,anddiscontinuousGalerkin)andmostofthemcaninprinciplebe directlyusedtosimulatetheMHDsystems,themainchallengeinnumericallysolvingthe MHDequationsisrelatedtothedivergence-freeconditiononthemagnetic.First,we notethattheMHDsystem(2.1)alongwith(2.3)isalreadyaclosedsetofeightevolution equations.Second,wenotethat r B =0isan involution insteadofaconstraint(see page119{128of[24]),becauseif r B =0isinitially( t =0),thensystem(2.1) guaranteesthat r B =0isforallfuturetime( t> 0).Unfortunately,most numericaldiscretizationsofMHDdonotpropagatesomediscreteversionof r B =0forward 21 intime.Ashasbeenshownrepeatedlyintheliterature,failuretoadequatelycontrolthe resultingdivergenceerrorscanleadtonumericalinstability(seee.g.,[6,47,56,60,71,75]). Toaddressthisissue,wewillmakeuseofthemagneticpotentialinthenumericalmethods describedinthiswork. Becauseitisdivergence-free,themagneticcanbewrittenasthecurlofamagnetic vectorpotential: B = r A : (2.15) Furthermore,wecanwritetheevolutionequationofthemagneticintheMHDsystems (2.1)incurlform: @ B @t + r ( B u )=0 ; (2.16) duetothefollowingrelation r ( u B B u )= r ( B u ) : (2.17) Usingthemagneticvectorpotential(2.15),evolutionequation(2.16)canbewrittenas r ˆ @ A @t +( r A ) u ˙ =0 : (2.18) Therelation(2.18)impliestheexistenceofascalarfunction suchthat @ A @t +( r A ) u = : (2.19) 22 Inordertouniquely(atleastuptoadditiveconstants)determinetheadditionalscalar function ,wemustprescribesomegaugecondition. Afterinvestigatingseveralgaugeconditions,Helzeletal.[34]foundthatonecanobtain stablesolutionsbyintroducingtheWeylgauge,i.e.,setting 0.Withthisgaugechoice, theevolutionequationforthevectorpotentialbecomes @ A @t +( r A ) u =0 ; (2.20) whichcanberewrittenasanon-conservativequasilinearsystem, @ A @t + N 1 @ A @x + N 2 @ A @y + N 3 @ A @z =0 ; (2.21) where N 1 = 2 6 6 6 6 6 4 0 u y u z 0 u x 0 00 u x 3 7 7 7 7 7 5 ; N 2 = 2 6 6 6 6 6 4 u y 00 u x 0 u z 00 u y 3 7 7 7 7 7 5 ; N 3 = 2 6 6 6 6 6 4 u z 00 0 u z 0 u x u y 0 3 7 7 7 7 7 5 : (2.22) Oneywithsystem(2.21){(2.22)isthatitisonlyweaklyhyperbolic[34].Inorder toseethisweakhyperbolicity,westartwiththeJacobianmatrixinsomearbitrary direction n =( n x ;n y ;n z ): n x N 1 + n y N 2 + n z N 3 = 2 6 6 6 6 6 4 n y u y + n z u z n x u y n x u z n y u x n x u x + n z u z n y u z n z u x n z u y n x u x + n y u y 3 7 7 7 7 7 5 : (2.23) 23 Theeigenvaluesofmatrix(2.23)are 1 =0 ; 2 = 3 = n u ; (2.24) andthematrixofrighteigenvectorscanbewrittenas R = 2 6 6 6 6 6 4 r (1) r (2) r (3) 3 7 7 7 7 7 5 = 2 6 6 6 6 6 4 n x n y u z n z u y u x ( u n ) n x k u k 2 n y n z u x n x u z u y ( u n ) n y k u k 2 n z n x u y n y u x u z ( u n ) n z k u k 2 3 7 7 7 7 7 5 : (2.25) Ifweassumethat k u k6 =0and k n k =1,thedeterminantofmatrix R is det( R )= u k 3 cos( )sin( ) ; (2.26) where istheanglebetween n and u .Inparticular,thereexistfourdegeneratedirections, =0 ;ˇ= 2 ;ˇ; and3 ˇ= 2,inwhichtheeigenvectorsareincomplete.Hence,thesystem(2.21) isonlyweaklyhyperbolic. 2.4Magneticpotentialin2D AspecialcaseofthesituationdescribedaboveistheMHDsystemin2D.Inparticular,what wemeanby2Disthatalleightconservedvariables, q =( ˆ;ˆ u ; E ; B ),canbenon-zero,but eachdependsononlythreeindependentvariables: t , x ,and y .Fromthepoint-of-viewof themagneticpotential,the2Dcaseismuchsimplerthanthefull3Dcase,duetothefact 24 thatthedivergence-freeconditionto r B = @B x @x + @B y @y =0 : (2.27) Itcanbereadilyseenthatsolving B 3 byanynumericalschemewillnothaveanyimpacton thesatisfactionofthedivergence-freecondition(2.27).Inotherwords,usingthemagnetic potentialtone B 3 isunnecessary.Instead,inthe2Dcase,wecanwrite B x = @A z @y and B y = @A z @x ; (2.28) whichinvolvesonlythethirdcomponentofthemagneticpotential,therebyelyre- ducingthemagnetic vector potentialtoa scalar potential.Consequently,theCTmethodin 2Dcanbetosolvinganadvectionequationforthethirdcomponentofthevector potential: @A z @t + u x @A z @x + u y @A z @y =0 : (2.29) Thishastheaddedbthat(2.29)is strongly hyperbolic,unlikeitscounterpartinthe 3Dcase. 2.5GeneralframeworkofCTalgorithm IntheunstaggeredCTmethodfortheidealMHDequations[35],Helzel etal. coupleda conservativevolumehyperbolicsolverfortheMHDequationswithanon-conservative volumesolverforthevectorpotentialequationtosolvetheidealMHDsystems.Using 25 theirapproachasthebasicframeworkforour2Dand3Dschemes,theworkinthispaper ismainlyfocusedonextendingthe3rd-ordervolumeCTmethodtoahigh-orderand computationallytCTmethodbutstillkeepingalltheadvantagesof theunstaggeredCTmethod.Inthissection,wewillsummarizethegeneralunstaggeredCT method. Assumethesemi-discreteformofMHDequations(2.1)hasageneralform Q 0 MHD ( t )= L ( Q MHD ( t ))(2.30) andthesemi-discreteformoftheevolutionequationofthemagneticpotential((2.21)for 3Dcaseand(2.29)for2Dcase)hasaform Q 0 A ( t )= H ( Q A ( t ) ; u ( t ))(2.31) wheretheabstractAisintroducedtodenotethevectorpotential A in3Dcaseorthe scalarpotential A z in2Dcase, Q MHD ( t )representsthegridfunctionattime t ofconserved quantitiesintheidealMHDsystemand Q A ( t )representsthegridfunctionofthemagnetic potential A or A z . Forsimplicity,weonlypresenttheCTschemescoupledwithforwardEulersteppingto solvetheMHDsystems.Whentheproblemissolvedfromthecurrentstate t = t n toanew state t = t n +1 ,asingletimestepoftheCTmethodconsistsofthefollowingsubsteps: 0. Startwith Q n MHD and Q n A (thesolutionsat t n ortheinitialconditionat t 0 ) 1. Buildtherighthandsidesofbothsemi-discretesystems(2.30)and(2.31)bysome 26 spatialdiscretizationsandupdatethesystemsby Q MHD = Q n MHD + t L ( Q n MHD ) Q n +1 A = Q n A + t H ( Q n A ; u n ) (2.32) where Q MHD =( ˆ n +1 ;ˆ u n +1 ; E ; B ), B meansthepredictedmagneticwithout satisfyingthedivergencefreeconstraintand E willbeupdatedbasedontheoption inStep4. 2. Correct B bythemagneticpotential Q n +1 A atnewstagebyadiscretecurloperator B n +1 = r Q n +1 A (2.33) 3. Setthetotalenergydensity E n +1 basedonthefollowingoptions: Option1: Keepthetotalenergyconserved: E n +1 = E (2.34) Option2: Keepthepressurethesameafterthecorrectionofmagnetic E n +1 = E + 1 2 ( k B n +1 k 2 k B k 2 )(2.35) Thishelpspreservethepositivityofthepressureforlowpressureproblems,which increasethestabilityofthenumericalsolver,althoughitenergyconser- vation. 27 Chapter3 FD-WENOwithconstrained transportforidealMHD Inthischapter,wedevelopaclassofhigh-orderweightedessentiallynon- oscillatoryschemesforsolvingtheidealMHDequationsin2Dand3D.Thephilosophyof thisworkistousethigh-orderWENOspatialdiscretizationswithhigh-orderSSP- RKtime-steppingschemes.Numericalresultsshowthatwithsuchmethodsweareableto resolvesolutionstructuresthatareonlyvisibleatmuchhighergridresolutionswithlower- orderschemes.ThekeychallengeinapplyingsuchmethodstoidealMHDistocontrol divergenceerrorsinthemagneticWeachievethisbyaugmentingthebasescheme withanovelhigh-orderconstrainedtransportapproachthatupdatesthemagneticvector potential.Thepredictedmagneticfromthebaseschemeisreplacedbyadivergence-free magneticthatisobtainedfromthecurlofthismagneticpotential.Thenon-conservative weaklyhyperbolicsystemthatthemagneticvectorpotentialissolvedusingaversion ofFD-WENOdevelopedforHamilton{Jacobiequations.Theresultingnumericalmethod isendowedwithseveralimportantproperties:(1)allquantities,includingallcomponents ofthemagneticandmagneticpotential,aretreatedaspointvaluesonthesamemesh (i.e.,thereisnomeshstaggering);(2)boththespatialandtemporalordersofaccuracyare Thischapterpreviouslyappearedas[22]:A.J.Christlieb,J.A.RossmanithandQ.Tang.Finitedierence weightedessentiallynon-oscillatoryschemeswithconstrainedtransportforidealmagnetohydrodynamics. J. Comput.Phys., 268:302{325,2014. 28 fourth-order;(3)nospatialintegrationormultidimensionalreconstructionsareneededin anystep;and(4)speciallimitersinthemagneticvectorpotentialupdateareusedtocontrol unphysicaloscillationsinthemagnetic Theoutlineofthischapterisasfollows.InSections3.1-3.4,wedescribetpartsof thenumericaldiscretizationoftheproposedmethod.InSection3.1wedetailthe5th-order FD-WENOspatialdiscretizationfortheMHDsystem.InSection3.2weoutlinethetime steppingtechniques.InSection3.3wedescribethespatialdiscretizationofthe2Dscalar potentialevolutionequationsandthe3Dvectorpotentialequations.Thenumericalcurl operatoranditspropertiesarethendiscussedinSection3.4,whichcompletesallthesteps ofourCTmethods.Theresulting2Dand3Dschemesareimplementedandtestedonseveral numericalexamplesinSection4.4. 3.1SpatialdiscretizationofidealMHD Inthissection,wedescribethesemi-discreteceweightedessentiallynon- oscillatoryschemethatcomprisesthe basescheme intheconstrainedtransportframework describedinSection2.5.OurmethodofchoiceistheFD-WENOmethoddevelopwedby JiangandWu[40].WewillrefertothismethodastheWENO-HCL 1 scheme.Inwhatfol- lows,wedescribethebasicWENO-HCLschemeinonespacedimensionfortheidealMHD equations.Wethendiscussthestraightforwardextensiontohigherdimensions. WewritetheMHDsystem(2.1)in1Dasfollows: @q @t + @ f ( q ) @x =0 ; (3.1) 1 WENO-HCL:=WeightedEssentiallyNon-OscillatoryforHyperbolicConservationLaws. 29 where q =( ˆ;ˆu x ;ˆu y ;ˆu z ; E ;B x ;B y ;B z ) T ; (3.2) f ( q )= ˆu x ;ˆu x u x + p + 1 2 k B k 2 B x B x ;ˆu x u y B x B y ;ˆu x u z B x B z ; u x E + p + 1 2 k B k 2 B x ( u B ) ; 0 ;u x B y u y B x ;u x B z u z B x T : (3.3) Forconvenience,wealsointroduce w =( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z ) T (3.4) todenotethevectorofprimitivevariables. DuetothehyperbolicityoftheMHDsystems,theJacobianmatrix @ f @q hasaspectral decompositionoftheform @ f @q = ; (3.5) where isthediagonalmatrixofrealeigenvalues, R isthematrixofrighteigenvectorsand L = R 1 isthematrixoflefteigenvectors. Weconsidertheproblemonauniformgridwith N +1gridpointsasfollows: a = x 1 2 > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > : 0if0 : 00 x 0 : 25 ; ( x 0 : 25) = 0 : 075if0 : 25 x 0 : 40 ; 2if0 : 40 x 0 : 60 ; (0 : 75 x ) = 0 : 075if0 : 60 x 0 : 75 ; 0if0 : 75 x 1 : 00 : (3.28) ThisproblemwasconsideredbyRossmanith[60],whereitwasalsousedtotestalimiter especiallydesignedtocontroloscillationsinthederivativeof q . ThesolutionsandtheirnumericalderivativescomputedbyWENO-HCLandWENO-HJ schemesarepresentedinFigure3.1.BothapproachesuseWENOreconstruction inthespatialdiscretizationandSSP-RK4forthetimeintegrator.WeuseaCFLnumberof 1.0andcomputethesolutionto t =4.Showninthisare3.1(a)thesolutionobtained bytheWENO-HCLschemeonameshwith N =300,3.1(b)thederivativeofthissolutionas computedwithafourth-ordercentralapproximation,3.1(c)thesolutionobtained bytheWENO-HJschemeonameshwith N =300,and3.1(d)thederivativeofthissolution ascomputedwithafourth-ordercentralapproximation.Althoughbothsolutions 40 (a) (b) (c) (d) Figure3.1:1Dadvectionequation.Showninthesepanelsare(a)thesolutionobtained bytheWENO-HCLschemewith(b)itsderivative,and(c)thesolutionobtainedbythe WENO-HJschemewith(d)itsderivative. Forinterpretationofthereferencestocolorin thisandallotheres,thereaderisreferredtotheelectronicversionofthisdissertation. agreewiththeexactsolutionverywell,thecomputedderivative @q @x (4 ;x )ofWENO-HCLis muchmoreoscillatorythanthatofWENO-HJ.TheproposedWENO-HJschemeisableto controlunphysicaloscillationsbothinthesolutionanditsderivative.Theresultofournew approachiscomparabletothatofexistingvolumeapproaches[34,35,60]. 41 3.3.22Dmagneticpotentialequation IntheCTframeworkdescribedinSection2.5,duringStep1wemustupdatethemagnetic potentialbysolvingadiscreteversionof @A z @t + u x ( x;y ) @A z @x + u y ( x;y ) @A z @y =0 ; (3.29) where,asdescribedinSection2.5,thevelocitycomponentsaregivenfunctionsfromthe previoustimestep(ortimestageinthecaseofhigher-ordertime-stepping).Because u x and u y aregiven,wecanview(3.29)asaHamilton{Jacobiequation: @A z @t + H x;y; @A z @x ; @A z @y =0 ; (3.30) withHamiltonian: H x;y; @A z @x ; @A z @y = u x ( x;y ) @A z @x + u y ( x;y ) @A z @y : (3.31) Tosolvethisequationwecandirectlyapplyatwo-dimensionalversionoftheWENO-HJ schemedescribedabove(justaswithWENO-HCL,the2DversionofWENO-HJissimply adirection-by-directionversionofthe1Dscheme).The2Dsemi-discreteWENO-HJcanbe writtenas dA z ij ( t ) dt = ^ H A z xij ;A z + xij ;A z yij ;A z + yij = u x ij A z xij + A z + xij 2 ! u y ij A z yij + A z + yij 2 ! + 1 A z + xij A z xij 2 ! + 2 A z + yij A z yij 2 ! ; (3.32) 42 where 1 =max i;j u x ij and 2 =max i;j u y ij : Theapproximations A 3 xij and A 3 yij arecalculatedwithformulasanalogousto(3.21)and (3.22). 1 and 2 arechosenasthemaximalvalueoverallgridpointsbasedonasimilar ideaoftheLax-FriedrichssplittinginSection3.1.Weremarkherethescheme(3.32) withthisglobal m canbetoodissipativeforcertainpurelinearHamilton{Jacobiequations. Anotherobviouschoiceistoevaluate m bytakingthemaximalvalueonthelocalstencil. Althoughthislocalversionof m canbemuchlessdissipativeforacertainpureHamilton{ Jacobiequation,innumericalexperimentsforMHDwethatthebetween thelocalandglobalapproachesarenegligible.Therefore,wewillonlypresentthenumerical resultsbytheglobalversionof m inthenumericalexamplessection(Section4.4). Exceptforthelastremainingdetailabouthowthediscretecurlof A z iscomputed(see Section3.4),thisversionoftheWENO-HJschemecoupledwiththeWENO-HCLscheme asthebaseschemecompletesour2DdiWENOconstrainedtransportmethod. Inthenumericalexamplessection,Section4.4,wewillrefertothisfullschemeasWENO- CT2D. 3.3.33Dmagneticpotentialequation Theevolutionequationofthemagneticpotential(2.21)in3Distlytfrom theevolutionequationofthescalarpotential(2.29)in2D,andhencetheschemediscussed inSection3.3.2cannotbeuseddirectly.AspointedoutinSection2.3,akeycultyisthe weakhyperbolicityofsystem(2.21).Helzel etal. [34]foundthattheweakhyperbolicityof 43 (2.21)isonlyanartifactoffreezingthevelocityintime,i.e.,thefullMHDsystemisstill stronglyhyperbolic.Furthermore,theyfoundthatthemagneticvectorpotentialsystem canbesolvedbyanoperatorsplitvolumeschemewithanadditionallimitingstrategy addedincertaindirections.Helzel etal. [35]handledtheweaklyhyperblicsystem(2.21) throughapath-conservingvolumeWENOschemewithoutappealtooperatorsplitting. Inordertoexplaintheschemeadvocatedinthiswork,wetakeinspirationfromthe operatorsplitmethodofHelzel etal. [34]andseparatesystem(2.21)intotwosub-problems: Sub-problem1: @A x @t + u y @A x @y + u z @A x @z =0 ; @A y @t + u x @A y @x + u z @A y @z =0 ; @A z @t + u x @A z @x + u y @A z @y =0 : (3.33) Sub-problem2: @A x @t u y @A y @x u z @A z @x =0 ; @A y @t u x @A x @y u z @A z @y =0 ; @A z @t u x @A x @z u y @A y @z =0 : (3.34) Weemphasisherethatourschemewillnotcontainanyoperatorsplitting,andthat thedivisionofthemagneticpotentialevolutionequationintotheabovetwosub-problems isonlyforthepurposeofexposition. Thesub-problemisacombinationofthreeindependentevolutionequations,eachof whichhasthesamemathematicalformasthe2Dscalarevolutionequation(3.29).Further- 44 more,thissub-problemisstronglyhyperbolic.Thus,atleastforthissub-problem,wecan simplyusetheWENO-HJschemedescribedinSection3.3.2tosolvethesethreeequations independently. Thesecondsub-problem,(3.34),isonlyweaklyhyperbolic.Forthisproblemweapplya WENOncediscretizationusingarithmeticaveragestosolutionderivatives atgridpoints.Forinstance,forthecomponent A x in(3.34),thesemi-discreteform becomes d dt A x ijk ( t )= u y ijk 0 @ A y xijk + A y + xijk 2 1 A + u z ijk A z xijk + A z + xijk 2 ! ; (3.35) where A xijk againusestheWENOreconstructiongivenby(3.21){(3.22).Thesemi-discrete formsfortheothercomponentsin(3.34)aresimilar. Notethatsemi-discreteformula(3.35)lacksthenumericaldissipationtermsfoundin (3.32).Thisisduetothefactthatsystem(3.34)(byitself)doesnotrepresentatransport equation.Therefore,theabovedescribeddiscretizationsfor(3.33)and(3.34)generallydo notintroducetnumericalresistivityinordertocontrolunphysicaloscillationsinthe magneticfora3Dproblem.Tobemoreprecise,whensolvingsystem(3.33), resistivityisintroducedfromtheWENOupwindingprocedure(seeformula(3.32)),but onlyin2ofthe3coordinatedirections(e.g.,for A z thereisresistivityinthe x and y -directions).Whensolvingsystem(3.34),noadditionalresistivityisintroduced (e.g.,see(3.35)).ThislackofnumericaldissipationwasalsopointedoutbyHelzel etal. [34,35];theyintroducedexplicitresistivitytermsintothemagneticvectorpotential equation.Wefollowasimilarapproachbymodifyingsub-problem(3.34)asfollows: 45 Sub-problem2withresistivity: @A x @t u y @A y @x u z @A z @x = " 1 @ 2 A x @x 2 ; @A y @t u x @A x @y u z @A z @y = " 2 @ 2 A y @y 2 ; @A z @t u x @A x @z u y @A y @z = " 3 @ 2 A z @z 2 : (3.36) Theseadditionaltermsgiveusresistivityinthemissingdirections(e.g.,theequation for A z nowhasanlresistivityterminthe z -direction).Intheaboveexpression,the resistivityistaketobeofthefollowingform: " 1 =2 1 x 2 + t ; (3.37) where0 ˝ 1issmallparameterthatcanbesettoensurethat " 1 remainsbounded as t ! 0 + , 1 isthesmoothnessindicatorof A x ,and isaconstantusedtocontrolthe magnitudeoftheresistivity. Inallthesimulationspresentedinthiswork,wetake =0 t and x havethesame orderofmagnitudeforalltheproblemsconsideredinthiswork).Thesmoothnessindicator 1 iscomputedasfollows: 1 ijk = a a + a + 1 2 ; (3.38) where a = ˆ + xA x xijk 2 ˙ 2 and a + = ˆ + xA x + xijk 2 ˙ 2 ; (3.39) 46 and istakentobe10 8 inallofournumericalcomputations.Here a and a + areused toindicatethesmoothnessof A x ineachofthe and+WENOstencils,respectively. Theresistivityparameters 2 and 3 intheotherdirectionscanbecomputedin analogousways. Thesmoothnessindicator i isdesignedsuchthattaresistivityisintro- ducedtoavoidspuriousoscillationsinthederivativesof A i when A i isnon-smooth,and high-orderaccuracyoftheschemeismaintainedwhen A i issmooth.Forthecasewhen @A x @x issmooth: A x xijk A x + xijk = O x 5 )and 1 ijk = O x 5 ) : Inthiscasetheresistivitytermin(3.36)willbeof O x 6 ),whichwillnotdestroy thespatialaccuracyofthescheme.Forthecasewhen @A x @x isnon-smooth: A x xijk A x + xijk = O (1)and 1 ijk ˇ 1 2 ; whichindicatesthatnumericalresistivityshouldbeadded.Forboththesmoothandnon- smoothcases,wenotethat 1 < 1 2 ,whichmeansthatforforwardEulertimestepping, in therangeof[0 ; 0 : 5]willguaranteethatthenumericalschemewillbestableuptoCFL=1. Forthefourth-order10-stageSSP-RK4time-steppingscheme,wefoundthat intherange of[0 : 02 ; 0 : 2]seemstosatisfactorilycontroltheunphysicaloscillationsin3Dproblems. Theconstrainedtransportmethodthatweadvocateinthisworkisamethodoflines approach,andthusisnotconsistentwiththeoperatorsplittingapproach.However,through numericalexperimentation,wediscoveredthatoperatorsplittingisnotnecessarytoobtain 47 accurateandstablesolutions,aslongastheaboveresistivitylimitingstrategyis includedinthetimeevolution.Inordertowriteoutthealschemeasadvocated,consider forbrevityonlytheequationinthemagneticvectorpotentialsystemwith resistivity: @A x @t u y @A y @x u z @A z @x + u y @A x @y + u z @A x @z = 1 @ 2 A x @x 2 : (3.40) Usingtheabovediscussionaboutsub-problems1and2asaguide,wearriveatthefollowing unsplitsemi-discreteformforthefull A x evolutionequation: dA x ijk ( t ) dt = u y ijk 0 @ A y xijk + A y + xijk 2 1 A + u z ijk A z xijk + A z + xijk 2 ! +2 1 A x i 1 jk 2 A x ijk + A x i +1 jk + t ! u y ijk A x yijk + A x + yijk 2 ! u z ijk A x zijk + A x + zijk 2 ! + 2 A x + yijk A x yijk 2 ! + 3 A x + zijk A x zijk 2 ! ; (3.41) where 2 =max i;j;k u y ijk ; and 3 =max i;j;k u z ijk : Thesemi-discreteformsfor A y and A z ofthesystemhaveanalogousforms.Forbrevitywe omittheseformulas. Exceptforthelastremainingdetailabouthowthediscretecurlof A iscomputed(see Section3.4),thisversionoftheWENO-HJschemecoupledwiththeWENO-HCLscheme 48 asthebaseschemecompletesour3DdiWENOconstrainedtransportmethod. Inthenumericalexamplessection,Section4.4,wewillrefertothisfullschemeasWENO- CT3D. Finally,wenotethattheresistivitytermsincludedinthesemi-discreteequation (3.41)with(3.37),(3.38),and(3.39)arespdesignedforsolvingtheidealMHD equations(2.1).Infutureworkwewillconsidernon-idealcorrections,includingphysical resistivityandtheHallterm.Inthesenon-idealcases,mowillhavetobemade totheresistivitytermsadvocatedinthiswork. 3.4Centraldiscretizationof r A DuringeachstageofourCTalgorithm,adiscretecurloperatorisappliedtothemagnetic potentialtogiveadivergence-freemagneticInthissectionwedescribetheapproach toapproximatethecurloperatoranddiscussitsimportantproperties. 3.4.1Curlin2D Welookforadiscreteversionofthe2Dcurlgivenby(2.28)ofthefollowingform: B x ij := D y ij A z and B y ij := D x ij A z ; (3.42) where D x and D y arediscreteversionsoftheoperators @ @x and @ @y .Inparticular,welook fordiscreteoperators D x and D y withthepropertythat r B ij := D x ij B x + D y ij B y = D x ij D y ij A z D y ij D x ij A z =0 ; (3.43) 49 whichmeansthatwesatisfyadiscretedivergence-freecondition.Inthesecondorderaccu- rate,unstaggered,CTmethodsdevelopedbyHelzel etal. [34],Rossmanith[60],andoth [75],theobviouschoicefor D x ij and D y ij aresecond-ordercentralInthis work,inordertomaintainhigh-orderaccuracy,wereplacesecond-ordercentral withfourth-ordercentral D x ij A := 1 x A i 2 j 8 A i 1 j +8 A i +1 j A i +2 j ; (3.44) D y ij A := 1 y A ij 2 8 A ij 1 +8 A ij +1 A ij +2 : (3.45) 3.4.2Curlin3D Welookforadiscreteversionofthe3Dcurlofthefollowingform: B x ijk := D y ijk A z D z ijk A y ; (3.46) B y ijk := D z ijk A x D x ijk A z ; (3.47) B z ijk := D x ijk A y D y ijk A x : (3.48) where D x , D y ,and D z arediscreteversionsoftheoperators @ @x , @ @y ,and @ @z .Inparticular, welookfordiscreteoperators D x , D y ,and D z withthepropertythat r B ijk := D x ijk B x + D y ijk B y + D z ijk B z = D x ijk D y ijk A z D x ijk D z ijk A y + D y ijk D z ijk A x D y ijk D x ijk A z + D z ijk D x ijk A y D z ijk D y ijk A x =0 ; (3.49) 50 whichmeansthatwesatisfyadiscretedivergence-freecondition.Toachievehigher-order accuracyweagainusefourth-ordercentral D x ijk A := 1 x A i 2 jk 8 A i 1 jk +8 A i +1 jk A i +2 jk ; (3.50) D y ijk A := 1 y A ij 2 k 8 A ij 1 k +8 A ij +1 k A ij +2 k ; (3.51) D z ijk A := 1 z A ijk 2 8 A ijk 1 +8 A ijk +1 A ijk +2 : (3.52) 3.4.3Importantproperties Forsmoothsolutions,thespatialaccuracyofouroverallschemewillbefourth-orderac- curate.ThisfactisvianumericalexperimentsinSection4.4.Furthermore,for solutionswithdiscontinuitiesinthemagneticthefourth-ordercentraldiscretizationof themagneticpotentialcurlwillintroducespuriousoscillations.However,aswedemonstrated vianumericalexperiments,weareabletocontrolanyunphysicaloscillationsinthemag- neticthroughthelimitingstrategythatwasdesignedinSection3.3fortheWENO-HJ scheme. Finally,wepointoutthefollowingpropertyoftheproposedscheme: Claim3.4.1. Theconstrainedtransportmethodasdescribedinthisworkgloballyconserves themagneticfromoneRunge{Kuttastagetothenext. Proof. Usingthesameideaastheproofoftheconservationof B inRossmanith[60],we canshowthetotalamountofeachcomponentof B canbemoonlythroughlossor gainacrossitsphysicalboundary.Thusthecomponentsofthemagneticareglobally conserved.Weomitthedetailsofproofhere. 51 3.5Numericalresults Inthissection,the2Dand3DWENO-CTschemesareappliedtoseveralMHDproblems. First,boththe2Dand3Dschemesaretestedonthe2Dand3DsmoothAlfvenwave problems,respectively.Theseproblemsareusedtodemonstratethattheproposedmethods arefourth-orderaccurate.Theschemeisalsotestedonarotatedshocktubeproblemin ordertoexaminetheshock-capturingabilityofthemethod,aswellastodemonstratethe successincontrollingdivergenceerrors.Alsoconsideredarethe2DOrszag-Tangvortex problemandthe2D,2.5D,and3Dversionsofthecloud-shockinteractionproblem. Foralltheexamplescomputedbelow,thegasconstantis =5 = 3andtheCFLnumber is3 : 0.Inthissection,wemakeexclusiveuseof Option1 andthusconservethetotalenergy. 3.5.1SmoothAlfvenwaveproblem Weconsider2Dand3DversionsofthesmoothAlfvenwaveproblemtodemonstrate theorderofaccuracyoftheproposedschemes. 3.5.1.12Dproblem Weperformaconvergencestudyofthe2Dschemeforthe2DsmoothAlfvenwaveproblem. Theinitialconditionsandthecomputationaldomainforthisversionaredescribedindetail inseveralpapers(e.g.,Section6.1.1onpage3818ofHelzel etal. [34]).The L 2 -errorsand L 1 -errorsofthemagneticandthemagneticscalarpotentialareshowninTables3.1. Fourth-orderconvergenceratesofallthequantitiesareobservedwhenCFL=3 : 0,which thetemporalandspatialorderofaccuracy. 52 Table3.1:Convergencestudyofthe2DAlfvenwaveforCFL=3 : 0.Shownarethe L 2 -errors and L 1 -errorsattime t =1ofthemagneticandmagneticpotentialascomputedby theWENO-CT2Dschemeatvariousgridresolutions. Mesh B x B y L 2 -error order L 1 -error order L 2 -error order L 1 -error order 16 32 9.73e-05 - 2.70e-04 - 2.13e-04 - 5.79e-04 - 32 64 4.07e-06 4.58 1.09e-05 4.64 9.34e-06 4.51 2.47e-05 4.55 64 128 2.02e-07 4.33 4.81e-07 4.50 4.62e-07 4.34 1.09e-06 4.50 128 256 1.17e-08 4.11 2.73e-08 4.14 2.60e-08 4.15 6.06e-08 4.17 256 512 7.15e-10 4.03 1.65e-09 4.05 1.55e-09 4.07 3.62e-09 4.07 512 1024 4.44e-11 4.01 1.02e-10 4.01 9.50e-11 4.03 2.22e-10 4.02 Mesh B z A z L 2 -error order L 1 -error order L 2 -error order L 1 -error order 16 32 2.83e-04 - 7.32e-04 - 3.03e-05 - 6.98e-05 - 32 64 9.39e-06 4.91 2.59e-05 4.82 1.37e-06 4.47 3.08e-06 4.50 64 128 2.88e-07 5.03 7.94e-07 5.03 7.07e-08 4.27 1.56e-07 4.30 128 256 9.30e-09 4.95 2.50e-08 4.99 4.14e-09 4.09 9.21e-09 4.09 256 512 3.40e-10 4.77 8.13e-10 4.94 2.54e-10 4.03 5.67e-10 4.02 512 1024 1.55e-11 4.46 3.61e-11 4.49 1.58e-11 4.01 3.53e-11 4.00 3.5.1.23Dproblem Wealsoperformaconvergencestudyofthe3Dschemeona3Dvariantofthesmooth Alfvenwaveproblem.Theinitialconditionsandthecomputationaldomainforthisversion aredescribedindetailinHelzel etal. [34](page3819inSection6.2.1).Theresultsofthe3D schemearepresentedinTables3.2.Fourth-orderaccuracyinallcomponentsare bythistestproblem. 53 Table3.2:Convergencestudyofthe3DAlfvenwaveforCFL=3 : 0.Shownarethe L 2 -errors and L 1 -errorsattime t =1ofallthemagneticandmagneticpotentialcomponentsas computedbytheWENO-CT3Dschemeonvariousgridresolutions. Mesh B x B y L 2 -error order L 1 -error order L 2 -error order L 1 -error order 16 32 32 6.92e-05 - 3.07e-04 - 1.16e-04 - 5.47e-04 - 32 64 64 2.73e-06 4.66 1.20e-05 4.68 4.98e-06 4.54 2.10e-05 4.70 64 128 128 1.28e-07 4.41 5.09e-07 4.56 2.45e-07 4.34 9.56e-07 4.46 Mesh B z A x L 2 -error order L 1 -error order L 2 -error order L 1 -error order 16 32 32 1.06e-04 - 5.47e-04 - 6.36e-06 - 3.04e-05 - 32 64 64 4.27e-06 4.64 1.80e-05 4.93 3.31e-07 4.27 1.28e-06 4.57 64 128 128 2.07e-07 4.37 8.21e-07 4.45 1.84e-08 4.17 6.72e-08 4.25 Mesh A y A z L 2 -error order L 1 -error order L 2 -error order L 1 -error order 16 32 32 1.42e-05 - 5.24e-05 - 1.53e-05 - 6.24e-05 - 32 64 64 5.82e-07 4.61 2.24e-06 4.55 6.31e-07 4.60 2.45e-06 4.67 64 128 128 2.88e-08 4.34 1.11e-07 4.33 3.06e-08 4.37 1.16e-07 4.41 3.5.22Drotatedshocktubeproblem Nextweconsidera1Dshocktubeproblemrotatedbyanangleof ina2Ddomain.The initialconditionsare ( ˆ;u ? ;u k ;u z ;p;B ? ;B k ;B z )= 8 > > > < > > > : (1 ; 0 : 4 ; 0 ; 0 ; 1 ; 0 : 75 ; 1 ; 0)if ˘< 0 ; (0 : 2 ; 0 : 4 ; 0 ; 0 ; 0 : 1 ; 0 : 75 ; 1 ; 0)if ˘> 0 ; (3.53) where ˘ = x cos + y sin and = x sin + y cos ; (3.54) and u ? and B ? arevectorcomponentsthatareperpendiculartotheshockinterface,and u k and B k arecomponentsthatareparalleltotheshockinterface.Theinitialmagnetic 54 (a) (b) Figure3.2:Therotatedshocktubeproblem.Showninthesepanelsare ˆ of(a)thebase WENOschemeand(b)theWENO-CT2Dscheme.30equallyspacedcontoursareshown foreachgraphintheranges ˆ 2 [0 : 1795 ; 1]. potentialis A z (0 ;˘ )= 8 > > > < > > > : 0 : 75 ˘ if ˘ 0 ; 0 : 75 + ˘ if ˘ 0 : (3.55) Wesolvethisproblemonthecomputationaldomain( x;y ) 2 [ 1 : 2 ; 1 : 2] [ 1 ; 1]witha 180 150mesh.Wetake =tan 1 (0 : 5).Zero-orderextrapolationboundaryconditionsare usedontheleftandrightboundaries.Onthetopandbottomboundariesalltheconserved quantitiesareextrapolatedinthedirectionparalleltotheshockinterface.Inaddition,tobe consistentwithzero-orderextrapolationboundaryconditionon B ,thelinearextrapolation ofthemagneticpotential A z isusedalongthecorrespondingdirections. ShowninFigure3.2arethedensitycontoursofthesolutionsascomputedusingthebase WENOschemeandtheWENO-CT2Dscheme.Fromthiswenotethatthesolution fromthebaseschemefromunphysicaloscillationsthatareduetotodivergence errorsinthemagneticwhiletheWENO-CT2Ddoesnotfromthisproblem.In 55 (a) (b) (c) (d) Figure3.3:Therotatedshocktubeproblem.Showninthesepanelsare(a)a1Dcutalong y =0of ˆ ascomputedwiththebaseWENOscheme,(b)azoomedinviewofthesame plot,(c)a1Dcutalong y =0of ˆ ascomputedwiththeWENO-CT2Dscheme,and(d)a zoomedinviewofthesameplot.Thesolidlineisahighlyresolved1Dsolution. Figures3.3and3.4wealsopresentacomparisonofthe2Dsolutionsalong y =0compared againsta1DWENO-HCLsolutiononameshwith N =5000.Fromthese itisagainclearthatthebaseWENOschemeproducesextraoscillationsandlargererrors, whilethesolutionbytheWENO-CT2Dschemeareingoodagreementwiththe1Dsolution. 56 (a) (b) (c) (d) Figure3.4:Therotatedshocktubeproblem.Showninthesepanelsare1Dcutsalong y =0 ofthemagnetic(a)perpendicular( B ? )and(b)parallel( B k )totheshockinterface ascomputedwiththebaseWENOschemeand1Dcutsalong y =0ofthemagnetic (c)perpendicular( B ? )and(d)parallel( B k )totheshockinterfaceascomputedwiththe WENO-CT2Ddscheme.Thesolidlineisahighlyresolved1Dsolution. 3.5.32DOrszag-Tangvortex NextweconsidertheOrszag-Tangvortexproblem,whichiswidelyconsideredasastandard testexampleforMHDintheliterature(e.g.[25,60,75,83]),sincethesolutionatlatetimes inthesimulationisquitesensitivetodivergenceerrors.Theproblemhasasmoothinitial 57 (a) (b) (c) (d) Figure3.5:TheOrszag-Tangvortexproblem.Showninthesepanelsare ˆ attimes(a) t =0 : 5,(b) t =2,(c) t =3,and(d) t =4ascomputedwiththeWENO-CT2Dschemeon a192 192mesh.15equallyspacedcontourswereusedforeachplot. conditiononthedouble-periodicbox[0 ; 2 ˇ ] [0 ; 2 ˇ ]: ˆ (0 ;x;y )= 2 ;u x (0 ;x;y )= sin( y ) ;u y (0 ;x;y )=sin( x ) ; (3.56) p (0 ;x;y )= ;B x (0 ;x;y )= sin( y ) ;B y (0 ;x;y )=sin(2 x ) ; (3.57) u z (0 ;x;y )= B z (0 ;x;y )=0 ; (3.58) 58 (a) (b) Figure3.6:TheOrszag-Tangvortexproblem.Showninthesepanelsare(a)thethermal pressureascomputedwiththeWENO-CT2Dschemeat t =3ona192 192meshand(b) asliceofthethermalpressureat y =0 : 625 ˇ . withtheinitialmagneticpotential: A z (0 ;x;y )=0 : 5cos(2 x )+cos( y ) : (3.59) Periodicboundaryconditionsareimposedonalltheboundaries.Astimeevolves,theso- lutionformsseveralshockwavesandavortexstructureinthemiddleofthecomputational domain. WesolvetheMHDequationsona192 192meshwiththeWENO-CT2Dscheme.In Figure3.5,wepresentthecontoursofdensityat t =0 : 5,2,3,and4.Asliceofthepressure at y =0 : 625 ˇ and t =3isshownintherightpanelofFigure3.6.Althoughtpapers displaythesolutionatttimesandresolutions,ourresultsareingoodagreement withthosegivenin[25,60,75,83].Wedidnotobservetoscillationsinanyofthe conservedquantities,evenwhenthesystemisevolvedouttolongtimes.Oursimulation wassuccessfullyrunto t =30withouttheintroductionofnegativepressureanywherein 59 thecomputationaldomain.Ontheotherhand,thebaseschemewithoutCTstepproduces negativepressuresataround t =4 : 0ona192 192mesh. 3.5.4Cloud-shockinteraction Finallyweconsidertheso-calledcloud-shockinteractionproblem,whichinvolvesastrong shockinteractingwithadensecloudthatisinhydrostaticequilibrium.Forthisproblem, weconsiderthe2D,2.5D,and3Dversionsoftheproposedmethod.The2Dand2.5D versionshavethesamephysicalsetup.However,2Dmeanstheproblemissolvedusingthe WENO-CT2Dscheme,and2.5Dmeansthatthemagneticandthepotentialaresolved usingtheWENO-CT3Dscheme,i.e.,withallthecomponentsofthemagneticupdated, althoughallthequantitiesarestillindependentof z . 3.5.4.12Dproblem InthisversionweconsideranMHDshockpropagatingtowardastationarybubble,withthe samesetupastheonein[60].Thecomputationdomainis( x;y ) 2 [0 ; 1] [0 ; 1]withw boundaryconditionat x =0andwboundaryconditionsonthethreeothersides.The initialconditionsconsistofashockinitializedat x =0 : 05: ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )(0 ;x;y ) = 8 > > > < > > > : (3 : 86859 ; 11 : 2536 ; 0 ; 0 ; 167 : 345 ; 0 ; 2 : 1826182 ; 2 : 1826182)if x< 0 : 05 ; (1 ; 0 ; 0 ; 0 ; 1 ; 0 ; 0 : 56418958 ; 0 : 56418958)if x> 0 : 05 ; (3.60) 60 (a) (b) Figure3.7:The2Dcloud-shockinteractionproblem.Showninthesepanelsareschlieren plotsattime t =0 : 06of(a)ln( ˆ )and(b)thenormof B .Thesolutionwasobtainedusing theWENO-CT2Dschemeona256 256mesh. andacircularcloudofdensity ˆ =10andradius r =0 : 15centeredat( x;y )=(0 : 25 ; 0 : 5). Theinitialscalarmagneticpotentialisgivenby A z (0 ;x;y )= 8 > > > < > > > : 2 : 1826182( x 0 : 05)if x 0 : 05 ; 0 : 56418958( x 0 : 05)if x 0 : 05 : (3.61) Thesolutioniscomputedona256 256mesh.ShowninFigure3.7areschlieren plotsofln( ˆ )and k B k .Ingeneral,thesolutionshowsgoodagreementwiththeonein [60],althoughtheWENO-CT2Dschemesshowssomehigher-resolutionfeatures.Thereis anoticeableextrastructurethatcanbeobservedaround x =0 : 75ofthedensityplot(see Figure3.7(a)).Wealsothatwhentheresolutionofthesolutionusingthesecond-order volumesolver[60]isdoubledtoameshof512 512,asimilarstructurestartsto appear 2 .AsimilarstructurecanbeobservedinthesolutionofDaiandWoodward(Figure 2 ThistestwasdoneusingfreelyavailableMHDCLAW[58]code. 61 18in[25])ona512 512mesh,althoughtheirproblemsetupisslightlyt(i.e.,they usedastationaryshockinsteadofastationarybubble).Fromanotherperspective,asimilar structurearound x =0 : 75canalwaysbeobservedintheschlierenplotsof k B k solvedby tmethods.Duetothosefacts,itisreasonabletobelievethatthesolutionshould consistofthisstructureandthatthehigh-ordersolverwithlessnumericaldissipationcan obtainthisstructurewithfewergridpointsthanlow-ordermethods. 3.5.4.22.5Dproblem Wealsoconsiderthe2DversionproblemwiththemagneticsolvedbyWENO-CT3D schemesoastocompareour3Dand2Dschemes.Wecallthisproblemthe2.5Dversiontobe consistentwith[34,35].Theproblemisinitializedin2.5Dwiththesameinitialconditionsas the2Dversion(3.60).However,asdiscussedinSection3.3,themagneticpotentialevolutions of2Dand3Daretlyt. In2.5D,sinceallquantitiesareindependentof z ,themagneticvectorpotential thefollowingevolutionequation @A x @t u y @A y @x u z @A z @x + u z @A x @y =0 ; @A y @t + u x @A y @x u x @A x @y u z @A z @y =0 ; @A z @t + u x @A z @x + u y @A z @y =0 : (3.62) Themagneticeld B x = @A z @y ;B y = @A z @x ; and B z = @A y @x @A x @y : (3.63) Forthisversionofthecloud-shockproblem,themagneticvectorpotentialisinitializedas 62 (a) (b) Figure3.8:The2Dcloud-shockinteractionproblem.Showninthesepanelsare B z attime t =0 : 06ona256 256meshassolvedwith(a)theWENO-CT2Dschemeand(b)the2.5D versionoftheWENO-CT3Dscheme.25equallyspacedcontourswereusedineachofthese panels. follows: A (0 ;x;y )= 8 > > > < > > > : (0 ; 2 : 1826182( x 0 : 05) ; 2 : 1826182( x 0 : 05))if x 0 : 05 ; (0 ; 0 : 56418958( x 0 : 05) ; 0 : 56418958( x 0 : 05))if x 0 : 05 : (3.64) Thebetweenthe2Dand2.5Dschemesisessentiallythat B z isnotcorrected intheCTstepinthe2Dscheme,whileinthe2.5Dweupdate B z bypartialderivativesof A x and A y asdescribedabove.Itisalsoworthwhiletopointoutthat B z ;z isidenticallyzero inthiscase,sotheupdateof B z in(3.63)willnotthedivergenceerror.Intheend, thetwoapproachesproduceverysimilarresults.ShowninFigure3.8arecontourplotsof B z usingthetwotapproaches.FortheWENO-CT3Dschemetheelimiter describedinSection3.3.3wasusedwith =0 : 1. Althoughthesemethodscompute B z intheverytways,thetwosolutionsin Figure3.8areingoodagreement.Thisresultgivesusthattheproposed3D 63 schemeisabletosolvetheproblemevenwithstrongshocks.Inaddition,thereareno toscillationsobservedinthe2.5Dsolution.Theseresultsalsocomparewellwith theresultsofthesplittevolumeCTapproachof[34]andbytheunsplitMOL volumeCTapproachof[35]. 3.5.4.33Dproblem Thelastproblemweconsiderisafully3Dversionofthecloud-shockinteractionproblem. Theinitialconditionsconsistofashockinitializedat x =0 : 05: ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )(0 ;x;y;z )= 8 > > > < > > > : (3 : 86859 ; 11 : 2536 ; 0 ; 0 ; 167 : 345 ; 0 ; 2 : 1826182 ; 2 : 1826182)if x< 0 : 05 ; (1 ; 0 ; 0 ; 0 ; 1 ; 0 ; 0 : 56418958 ; 0 : 56418958)if x> 0 : 05 ; (3.65) andasphericalcloudofdensity ˆ =10andradius r =0 : 15centeredat( x;y;z )= (0 : 25 ; 0 : 5 ; 0 : 5).Theinitialconditionsforthemagneticpotentialarethesameas(3.64). Thesolutioniscomputedonthedomainof[0 ; 1] 3 .wboundaryconditionsareimposed at x =0andwboundaryconditionsareusedonallothersides.Thesolutioniscom- putedona128 128 128meshusingtheWENO-CT3Dscheme.InFigure3.9weshowthe densityofthesolutionat t =0 : 06,whichisingoodagreementwiththesolutionin[34,35], althoughoursolutionshowslessoscillationsandhigher-resolutioncomparedwithprevious work.Weagainobserveanextrastructureinthedensityplot,sameasinoursolutionofthe 2Dproblem.HeretheelimiterasdescribedinSection3.3.3wasusedwith =0 : 1. 64 Figure3.9:The3Dcloud-shockinteractionproblem.Showninthisplotisln( ˆ )ascomputed withtheWENO-CT3Dschemeona128 128 128mesh. 65 Chapter4 Positivity-preservinglimiterforideal MHD Inthischapter,weutilizethemaximum-principle-preservinglimitingtechnique,origi- nallydesignedforhighorderWENO-HCLmethodsforscalarhyperbolicconservationlaws, todevelopapositivity-preservinglimiterfortheFD-WENOmethodsfortheidealMHD equations.Theresultingscheme,undertheconstrainedtransportframework,canachieve highorderaccuracy,adiscretedivergence-freeconditionandpositivityofthenumerical solutionsimultaneously. Theoutlineofthischapterisasfollows.InSection4.1and4.2,wepresentpositivity- preservinglimitersforthe1Dandmulti-DMHDequations.Inaddition,wealsodetailthe temporaldiscretizationinthisworkinSection4.2.Theproposedschemesareimplemented andtestedonseveral1D,2Dand3DnumericalexamplesinSection4.4. 4.11Dcase Inthissection,wedescribeourpositivity-preservinglimiterona1DMHDsystem.The divergence-freecondition r B =0in1Dcaseisequivalentto B x =constant.Sincethe Thischapteristoappearin[21]:A.J.Christlieb,Y.Liu,Q.TangandZ.Xu.Positivity-preserving weightedENOschemeswithconstrainedtransportforidealmagnetohydrodynamicequations. SIAMJ.Sci.Comput., toappear,2015. 66 WENOhyperbolicconservationlawsolver(WENO-HCL)in[39,40]withoutCTapproaches willproduceasolutionwithconstant B x ,weuseitasourMHDbaseschemein1D,towhich weapplyapositivity-preservinglimiter. TheMHDequationsin1D(3.1)hasbeendiscussedinSection3.1.Again,welet q i ( t ) bethenumericalsolutionatthegridpoint x i .TheWENO-HCLschemes solve(3.1)byaconservativeform: dq i ( t ) dt + 1 x ^ F i + 1 2 ^ F i 1 2 =0 ; (4.1) where ^ F i + 1 2 isasahigh-ordernumericalconstructedbyWENO-HCL.The designof ^ F i + 1 2 hasbeenaddressedinSection3.1.Onenumericalyisthatthewave speedsoftheMHDsysteminvolvetheterm 1 ˆ .Toavoidthepossibilityofaniwave speedduringthecomputation,weassumethereisasmalllowerbound 0 forbothdensity andpressureintheexactsolutionoftheproblemweconsidered. Thesemi-discretizedequation(4.1)canbefurtherdiscretizedintimebyhigh-ordertime integrators.WhileourproposedschemecanbeappliedwithanySSP-RKmethod,wetake thethird-orderSSP-RKmethod(1.13)asanillustrativeexample.Herewelet L ( q n i )= 1 x ( ^ F n i + 1 2 ^ F n i 1 2 ) : (4.2) Ifweuse ^ F n i + 1 2 , ^ F (1) i + 1 2 and ^ F (2) i + 1 2 todenotethenumericalreconstructedbasedon q n , q (1) and q (2) ,thestageofRKdiscretization(1.13)canberewrittenas q n +1 i = q n i ( ^ F rk i + 1 2 ^ F rk i 1 2 ) ; (4.3) 67 where = t x ; ^ F rk i + 1 2 = 1 6 ^ F n i + 1 2 +4 ^ F (2) i + 1 2 + ^ F (1) i + 1 2 : (4.4) ^ F rk i + 1 2 canbeviewedasalinearcombinationofhigh-ordernumericalfromt stages.Followingtheideasin[78,79],weneedtomodifythenumerical ^ F rk i + 1 2 bya positivity-preservingtodesignahigh-orderpositivity-preservingMHDscheme. Cheng etal. [19]provedthesimpleLax{Friedrichsnumericalcoupledwithforward Eulertimediscretizationispositivity-preservingforthe1DMHDequations(3.1)underthe restrictionCFL 0 : 5.WhentheLax{Friedrichsschemeisusedtosolvethehigh-order solution q n from t n to t n +1 ,wehave ^ q n +1 i = q n i ( ^ f i + 1 2 ^ f i 1 2 ) ; (4.5) where^ q n +1 i isintroducedtodenotethelow-ordersolutionat x i and t = t n +1 ,andthe Lax{Friedrichsisformulatedas ^ f i + 1 2 = 1 2 f ( q n j +1 )+ f ( q n i ) ( q n j +1 q n i ) ; (4.6) wherethemaximalwavespeed isby =max i j u x j + c f : (4.7) Here c f isthefastspeedoftheMHDsystem,seeSection2.2. 68 Thedensityandpressurecomputedbythescheme(4.5)satisfy 8 > > > < > > > : ^ ˆ n +1 i > 0 ; ^ p n +1 i > 0 : (4.8) Wecanusethesolution^ q n +1 toethenumericallowerboundsforthedensity andpressureforthehigh-ordersolution q n +1 ,whichare n +1 ˆ =min i (^ ˆ n +1 i ; 0 ) ; (4.9) n +1 p =min i (^ p n +1 i ; 0 ) : (4.10) Throughoutthesimulationsforthiswork,wetake 0 =10 13 .Itcanbecertainlytakenas asmallernumberifitisrequiredbytheproblemandallowedbythemachineprecision. Following[78,79],toguaranteethepositivityofthehigh-ordersolutionsbytheWENO scheme(4.3),weneedtoamoofthenumericalasfollows: ~ F i + 1 2 = i + 1 2 ( ^ F rk i + 1 2 ^ f i + 1 2 )+ ^ f i + 1 2 ; (4.11) wherethelimitingparameter i + 1 2 2 [0 ; 1].Weseekacombinationof i + 1 2 ,suchthatthe solutionssatisfy 8 > > > < > > > : ˆ n +1 i n +1 ˆ ; p n +1 i n +1 p : (4.12) Ourpositivity-preservinglimitingtechniquefollowsatwo-stepprocedure.First,asout- 69 linedbelow,wedescribeastrategytoguaranteethecomputeddensitypositive.Tofacilitate thediscussion,wedenotetheofthedensityin ^ f as ^ f ˆ ,whereas f ˆ and ~ f ˆ arethecorrespondingcomponentsin ^ F rk and ~ F ,respectively. Topreservepositivedensity,weneedtoupperbounds ˆ 1 2 ;I i ofthelimitingpa- rameters i 1 2 ateachcell I i ,suchthat,foranycombination( i 1 2 ; i + 1 2 ) 2 [0 ; ˆ 1 2 ;I i ] [0 ; ˆ + 1 2 ;I i ],thefollowinginequalityholds: ˆ n +1 i ( i 1 2 ; i + 1 2 )= ˆ n i ( ~ f ˆ i + 1 2 ~ f ˆ i 1 2 ) n +1 ˆ ; (4.13) where ~ f ˆ i + 1 2 = i + 1 2 ( f ˆ i + 1 2 ^ f ˆ i + 1 2 )+ ^ f ˆ i + 1 2 ; Sinceweknow^ ˆ n +1 i = ˆ n i ( ^ f ˆ i + 1 2 ^ f ˆ i 1 2 ),(4.13) isequivalentto ^ ˆ n +1 i ( i + 1 2 ( f ˆ i + 1 2 ^ f ˆ i + 1 2 ) i 1 2 ( f ˆ i 1 2 ^ f ˆ i 1 2 )) n +1 ˆ : (4.14) Duetothepositivity-preservingpropertyoftheschemeandtheof n +1 ˆ (4.9),wehave^ ˆ n +1 i n +1 ˆ .Thus,theinequality(4.14)canberewrittenas, i 1 2 ( f ˆ i 1 2 ^ f ˆ i 1 2 ) i + 1 2 ( f ˆ i + 1 2 ^ f ˆ i + 1 2 ) n +1 ˆ ^ ˆ n +1 i (4.15) withtheright-handside n +1 ˆ ^ ˆ n +1 i 0.Forabbreviation,weintroduceanotation f i + 1 2 = f ˆ i + 1 2 ^ f ˆ i + 1 2 . Followingthesameideain[78,79],wewilldeterminetheupperboundsoftheparameter i 1 2 byacase-by-casediscussionbasedonthesignsof f i 1 2 and f i + 1 2 .Inparticular, wedecoupletheinequalities(4.15)basedonthefollowingfourcases: 70 If f i 1 2 0and f i + 1 2 0,then ˆ 1 2 ;I i ; ˆ + 1 2 ;I i )=(1 ; 1) : If f i 1 2 0and f i + 1 2 > 0,then ˆ 1 2 ;I i ; ˆ + 1 2 ;I i )= 0 @ 1 ; min 0 @ 1 ; n +1 ˆ ^ ˆ n +1 i f i + 1 2 1 A 1 A : If f i 1 2 < 0and f i + 1 2 0,then ˆ 1 2 ;I i ; ˆ + 1 2 ;I i )= 0 @ min 0 @ 1 ; n +1 ˆ ^ ˆ n +1 i f i 1 2 1 A ; 1 1 A : If f i 1 2 < 0and f i + 1 2 > 0, { iftheinequality(4.15)iswith( i 1 2 ; i + 1 2 )=(1 ; 1)then ˆ 1 2 ;I i ; ˆ + 1 2 ;I i )=(1 ; 1) : { otherwise,wechoose ˆ 1 2 ;I i ; ˆ + 1 2 ;I i )= 0 @ n +1 ˆ ^ ˆ n +1 i f i 1 2 f i + 1 2 ; n +1 ˆ ^ ˆ n +1 i f i 1 2 f i + 1 2 1 A : Thisprocedurehasbeendiscussedin[78,79].Itiseasytoshowwhen( i 1 2 ; i + 1 2 ) 2 [0 ; ˆ 1 2 ;I i ] [0 ; ˆ + 1 2 ;I i ]withthebounds ˆ 1 2 ;I i obtainedbytheabovestrategythat,the inequality(4.15)holds,i.e.,thedensity ˆ n +1 i ispositiveateachgrid x i .Wethisset 71 as S ˆ;I i : S ˆ;I i =[0 ; ˆ 1 2 ;I i ] [0 ; ˆ + 1 2 ;I i ] : (4.16) Wenextdescribeastrategytoobtainpositivepressure.Firstwediscusssomeproperties ofthepressurefunction p ( q )=( 1) E 1 2 ( ˆu x ) 2 +( ˆu y ) 2 +( ˆu z ) 2 ˆ 1 2 ( B x ) 2 +( B y ) 2 +( B z ) 2 : (4.17) Wenotethepressurefunctionisconcavewithrespectto q =( ˆ;ˆu x ;ˆu y ;ˆu z ; E ;B x ;B y ;B z ) T . Similarasthefunction ˆ n +1 i ( i 1 2 ; i + 1 2 ),wecanafunction p n +1 i ( i 1 2 ; i + 1 2 )as follows: p n +1 i ( i 1 2 ; i + 1 2 ):= p ( q n +1 i ( i 1 2 ; i + 1 2 )) : (4.18) Weneedthefollowinglemmatoconstructthelimiter. Lemma4.1.1. Thepressurefunction p q n +1 i ! 1 +(1 ) ! 2 p q n +1 i ! 1 +(1 ) p q n +1 i ! 2 (4.19) forany 2 [0 ; 1] and ! 1 ; ! 2 2 S ˆ;I i . Theproofofthislemmaisstraightforward,aslongasweusetheconcavepropertyof p ( q )andnotethatthesolution q n +1 i isalinearfunctionofitslimitingparameters,i.e., q n +1 i ! 1 +(1 ) ! 2 = q n +1 i ! 1 +(1 ) q n +1 i ! 2 : 72 AsimilarlemmafortheEulerequationshasbeenuseinthepast[20,78]. Wewanttoidentifyasubsetoftheset S ˆ;I i ,denotedby S p;I i ,suchthat p n +1 i ( i 1 2 ; i + 1 2 ) ispositive,i.e., S p;I i = f ( i 1 2 ; i + 1 2 ) 2 [0 ; ˆ 1 2 ;I i ] [0 ; ˆ + 1 2 ;I i ]: p n +1 i ( i 1 2 ; i + 1 2 ) n +1 p g : (4.20) DuetoLemma4.1.1, S p;I i isaconvexset.Todetermine S p;I i ,wecanonlyfocusonits vertices. Ifwedenotethefourverticesof S ˆ;I i tobe P k 1 ;k 2 =( k 1 ˆ 1 2 ;j ;k 2 ˆ + 1 2 ;j ) ; with k 1 , k 2 being0or1,similarlywecantheverticesof S p;I i tobe ^ P k 1 ;k 2 : Wedeterminethe ^ P k 1 ;k 2 basedonthefollowingstrategy.For( k 1 ;k 2 ) 6 =(0 ; 0),if p n +1 i ( P k 1 ;k 2 ) n +1 p ,welet ^ P k 1 ;k 2 = P k 1 ;k 2 ;otherwiseweascalarparameterof r suchthat p n +1 i ( rP k 1 ;k 2 ) n +1 p andlet ^ P k 1 ;k 2 = rP k 1 ;k 2 .Theresultingthreevertices ^ P k 1 ;k 2 withtheorigin(0 ; 0)form S p;I i . Next,wecanidentifyarectangleinside S p;I i denotedby R ˆ;p;I i =[0 ; 1 2 ;I i ] [0 ; + 1 2 ;I i ] ; (4.21) where 1 2 ;I i =min k 2 =0 ; 1 ( ^ P 1 ;k 2 ) ; + 1 2 ;I i =min k 1 =0 ; 1 ( ^ P k 1 ; 1 ) : (4.22) 73 Afterrepeatingthisprocedureforall j ,welet i + 1 2 = + 1 2 ;I i ; 1 2 ;I i +1 ) ; (4.23) andthisourdiscussionforthe1Dlimiter. Remark4.1.2. Thelimitingtechniquehereisusedonlytoguaranteethepositivityofthe solutionatthestageofRKmethods.Ifthereisnegativedensityorpressureinthe intermediatestage,inthisworkwetaketheabsolutevalueofthedensityandpressurein thecodewhereapositivesolutionisrequired.Theplacethatneedsapositivesolution istoestimatethespeedwavesofthesystem.Forinstance,thespeedofsoundistakenas c = p j p j = j ˆ j intheintermediatestage.Thesecondplacerequiringapositivesolutionis toestimatetheeigenvectorsoftheJacobianmatrixofthefunction.Thosetreatments willnotdegradetheorderofaccuracy,becausetheWENOalgorithmonlyneedsanestimate ofthelocaleigenvaluesandeigenvectorsandwealwaysusethetruesolutiontocompute thenumericalevenwhenitbecomesnegativeintheintermediatestage.Ontheother hand,wealsoremarkthatthelimitercanbeappliedtoeachstagewhenthepositivityinthe intermediatestageisrequired. Remark4.1.3. Fromthelimitingsteps,wecanseetheoverallschemehasaCFLconstraint of0.5,whichissameastheLax{Friedrichsscheme.Whentheabovelimitingtechniqueis appliedtotheintermediatestages,thereisnoextrarestrictionbecausethetimestepofthe intermediatestageistypicallynogreaterthan t . Remark4.1.4. Onenumericalistosatisfy p ( rP k 1 ;k 2 ) n +1 p .Thiscanbedone bysolvingaroot r fortheequation p ( rP k 1 ;k 2 )= n +1 p .Throughasimplederivation,itcan 74 beeasilyshownthatthesolution q n +1 i ( rP k 1 ;k 2 ) q n +1 i ( rP k 1 ;k 2 )= rq n +1 i ( P k 1 ;k 2 )+(1 r )^ q n +1 i ; where ^ q n +1 i isagainusedtodenotethesolutionsolvedbytheder ^ f i + 1 2 .This propertyisindependentfromthedimension,whichmakesitnaturallyextendableforthe multi-Dcases.Moreimportantly, q n +1 i ( P k 1 ;k 2 ) and ^ q n +1 i arebothcomputationallycheapto evaluate.Sowith q n +1 i ( P k 1 ;k 2 ) and ^ q n +1 i known,wecansolvearoot r fortheequation p ( rq n +1 i ( P k 1 ;k 2 )+(1 r )^ q n +1 i )= n +1 p : IntheMHDequationcase,thisequationisacubicfunctionof r ingeneral.Wenotethat thereexistsatleastonerootintheinterval [0 ; 1] ,whichcanalwaysbefoundbyNewton iteration.However,intheimplementation,weonlyusedasimplebisectionmethodwitha maximumof10iterationstotheroot,becauseourpurposeistoobtainapositivepressure p ( rq n +1 i ( P k 1 ;k 2 )) insteadofanaccurate r .Duringthenumericalsimulations,we foundthectofthenumberofiterationsonthesolutionqualityandaccuracyisnegligible. Asimilarapproachtoalimitingparametercanbefoundinthedesignofthepositivity- preservingMHDscheme[5]. 4.2Multi-Dcase Inthissection,wedescribethepositivity-preservinglimiterinthemulti-Dcase.To controlthedivergenceerror,ourbaseschemeistakenastheWENO-CTschemeoutlinedin Chapter3.Inthediscussionbelow,weonlyfocusonthe2Dcase,keepinginmindthatthe 75 extensionto3Dcaseisquitestraightforward. The2DMHDsystem(2.1)canberewrittenas: @q @t + @ f ( q ) @x + @ g ( q ) @y =0 : (4.24) Weneedtosolve(4.24)togettheupdatefortheconservedquantities.IfageneralRK methodisusedasthetimeintegrator,theWENO-HCLschemesolves(4.24)byaconservative form: q n +1 ij = q n ij x ( ^ F rk i + 1 2 j ^ F rk i 1 2 j ) y ( ^ G rk ij + 1 2 ^ G rk ij 1 2 ) ; (4.25) where ^ F rk and ^ G rk arelinearcombinationsofhigh-ordernumericalfromthreeRK stagesand = t .Let ^ f i + 1 2 j and ^ g ij + 1 2 againbethest-orderLax{Friedrichs Thenwemodifythehigh-ordernumerical ^ F rk and ^ G rk bytheLax{Friedrichs ^ f i + 1 2 j and ^ g ij + 1 2 toachievethepositivityofthesolution,i.e., ~ F i + 1 2 j = i + 1 2 j ( ^ F rk i + 1 2 j ^ f i + 1 2 j )+ ^ f i + 1 2 j ; (4.26) ~ G ij + 1 2 = ij + 1 2 ( ^ G rk ij + 1 2 ^ g ij + 1 2 )+ ^ g ij + 1 2 : (4.27) Foreachgrid x ij ,weuseatwo-stepstrategysimilartothe1Dcase.Firstly,wefollowa similarstrategytoguaranteethecomputeddensitypositive,i.e.,tofourupperbounds ˆ L;I ij , ˆ R;I ij , ˆ U;I ij and ˆ D;I ij suchthatforany( i 1 2 j ; i + 1 2 j ; ij 1 2 ; ij + 1 2 ) 2 S ˆ;I ij , 76 where S ˆ;I ij =[0 ; ˆ L;I ij ] [0 ; ˆ R;I ij ] [0 ; ˆ D;I ij ] [0 ; ˆ U;I ij ] ; (4.28) thecomputeddensity ˆ n +1 ij ( i 1 2 j ; i + 1 2 j ; ij 1 2 ; ij + 1 2 ) n +1 ˆ : (4.29) Secondly,wecanarectangularset R ˆ;p;I ij =[0 ; L;I ij ] [0 ; R;I ij ] [0 ; D;I ij ] [0 ; U;I ij ],whichisasubsetof S ˆ;I ij ,suchthatforany( i 1 2 j ; i + 1 2 j ; ij 1 2 ; ij + 1 2 ) 2 R ˆ;p;I ij ,wehave, p n +1 ij ( i 1 2 j ; i + 1 2 j ; ij 1 2 ; ij + 1 2 ) n +1 p : (4.30) Here n +1 ˆ and n +1 p arethe2Dlowerboundswithsimilarasthe1Dcase(4.9) and(4.10).Thestrategytotheset R ˆ;p;I ij issimilartotheEulerequationscase[78]. Weomitthedetailshere.Animportanttfromthe1DcaseinSection4.1is,thestep 4.22in1Dcasebecomes L;I ij =min k 2 ; 3 ; 4 =0 ; 1 ( ^ P 1 ;k 2 ;k 3 ;k 4 ) ; R;I ij =min k 1 ; 3 ; 4 =0 ; 1 ( ^ P k 1 ; 1 ;k 3 ;k 4 ) ; (4.31) D;I ij =min k 1 ; 2 ; 4 =0 ; 1 ( ^ P k 1 ;k 2 ; 1 ;k 4 ) ; U;I ij =min k 1 ; 2 ; 3 =0 ; 1 ( ^ P k 1 ;k 2 ;k 3 ; 1 ) : (4.32) 77 Afterrepeatingthisprocedureforallnodes( i;j ),welet i + 1 2 j = R;I ij ; L;I i +1 j ) ; (4.33) ij + 1 2 = U;I ij ; D;I ij +1 ) : (4.34) Thiswholeprocedurewillproduceanumericalsolutionwithpositivedensityandpressure afterStep1inCTframework.FollowedbyStep2and3with Option2 ,weachievehigh orderaccuracy,adiscretedivergence-freeconditionandpositivityofthenumericalsolution simultaneously.TheoverallschemesharesthesameCFLconstraintasthelow-orderLax- Fridrichsscheme.TherearenoextraCFLrestrictionsfromthelimitingprocess. Aspointedoutin[19],thereisstillnorigorousproofthattheLax{Friedrichsschemeor anyotherschemeispositivity-preservinginthemutli-Dcasewhenthedivergence- freeconstraintisconsidered.Inthiswork,westillusetheLax{Friedrichsscheme asthelow-ordercorrectionschemeforthemulti-Dcases.Like[19],wetakeCFL 0 : 5asthe constraintforthepositivity-preservingLax{Friedrichsschemeinthemulti-Dcases.Onthe otherhand,ourlimitingtechniqueisindependentfromthechoiceofthelow-orderscheme. Theoverallschemewillbeimprovedaslongasaprovablepositivity-preservingschemeis foundandusedasthebuildingblock. 4.3Temporaldiscretization Inthissectionwedescribesomedetailsoftime-steppingtechniquesusedinthiswork.In Chapter3,weuseSSP-RK4asourtimeintegratorandalsocorrectthemagneticeach timestep.Herewemakesomeadjustmentssothatthelimitingschemesbecomemore 78 t. First,wemakeexclusiveuseof Option2 inthiswork,soastopreservethepositivity ofthepressureafterthemagneticiscorrected.ThisisacommontechniqueintheCT frameworkforproblemsinvolvinglow- plasma[10,75].Underthisoption,ifthedensity andpressureafterStep1arepositive,theywillremainpositiveintheoverallcomputation. Therefore,withtheproposedlimiter,thenumericalsolutionsmaintainpositivethroughout thewholeCTsteps. Second,sincethelimiterisonlyappliedtothestageoftheRKmethodsinthis work,theRKmethodswithlargerCFLnumber(suchasSSP-RK4methodusedinChapter 3)isnotverystablebythistreatment.Soinstead,weusetraditionalSSP-RK3methodas thetestedintegratorwithaCFLnumberof0 : 5beingused. Anotheroftheschemesliesintheimplementationofthecorrectionsteps(Steps 2and3).Weproposetoperformthecorrectionstepsonlyattheendofeachtimestep t n insteadofeachstageofRKmethodsinChapter3.Thankstothismowecan simplyfocusonthestageofthesolutionwhenimplementingthelimitingtechnique inSection4.1.Numericalresultsshownegligiblebetweenthetwoapproaches whenSSP-RK3time-steppingisused.However,wenotethatthismomayresult inaccumulationofthedivergenceerrorespeciallyforRKmethodswithlargestagenumbers. Forthosetimesteppingschemes,thiskindofmoisnotrecommendedandthe correctionstepshavetobeperformedateachstage.Thus,ifSSP-RK4methodisusedfor thoselowdensityandpressureproblem,boththeproposedlimiterandthecorrectionsteps (Steps2and3)havetobeappliedateachtimestep. 79 4.4Numericalexamples Inthissection,weperformnumericalsimulationswiththeresultingpositivity-preserving schemesin1D,2Dand3D.SSP-RK3schemeservesasthetimeintegratorinalltheexamples, whereastheerWENO-HCLschemeisusedforsolvingthebaseMHD equationsintexamples.Inmulti-D,ourfourth-orderCTmethodisusedtoobtain adivergence-freemagneticUnlessotherwisestated,thegasconstantis =5 = 3and theCFLnumberis0 : 5. 4.4.1Testcasesin1D Inthissubsection,wetestourpositivity-preservingschemewithseveral1DMHDexamples. Wenotethatforallthecasespresentedinthissubsection,negativepressureordensityis observedifthebaseMHDschemeisappliedwithoutapositivity-preservinglimiter.Here, thebaseMHDschemeisaWENO-HCLscheme. 4.4.1.1Vacuumshocktubetest Weconsidera1Dvacuumshocktubeproblem.Thisexampleisusedtodemonstrate thatourpositivity-preservingMHDsolvercanhandleverylowdensityandpressure.The initialconditionis ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )= 8 > > > < > > > : (10 12 ; 0 ; 0 ; 0 ; 10 12 ; 0 ; 0 ; 0)if x< 0 ; (1 ; 0 ; 0 ; 0 ; 0 : 5 ; 0 ; 1 ; 0)if x> 0 : (4.35) Itissimilartothevacuumshocktubeproblemin[76].Thecomputationaldomainis [ 0 : 5 ; 0 : 5]andzero-orderextrapolationboundaryconditionsareused.ShowninFigure 80 (a) (b) Figure4.1:Thevacuumshocktubeproblem.Showninthesepanelsareplotsattime t =0 : 1 of(a) ˆ and(b)thethermalpressure.Thebluecircleisasolutionsolvedonameshwith N =200.Thesolidlineisahighlyresolvedsolution. 4.1arethedensityandpressureofthesolutiononameshwith N =200andthehighly resolvedsolutionwith N =2000.Wecanobservethesolutionsoflowresolutionandhigh resolutionareingoodagreement. 4.4.1.2TorsionalAlfvenwavepulse WealsoconsiderthetorsionalAlfvenwavepulseproblem[9,19].Theinitialconditionis ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )=(1 ; 10 ; 10cos ˚; 10sin ˚; 0 : 01 ; 10cos ˚; 10sin ˚; 0) ; (4.36) where ˚ = ˇ 8 (tanh( 0 : 25+ x +1))(tanh( 0 : 25 x +1))and =0 : 005.Thecomputational domainis[ 0 : 5 ; 0 : 5]andperiodicboundaryconditionsareused.Inthistestproblem,the initialpressureissosmallthattheproblemisverysensitivetothedissipationintroduced bynumericalschemes.Further,theexistenceofastrongtorsionalAlfvenwavediscontinuity makestheproblemulttosimulate.Inthesimulationwithoutthelimiter,thebase WENO-HCLintroducedanegativepressureinafewtimestepsandthesolutionsbecome unphysicalimmediately.Withthelimiter,ourschemecansimulatetheproblemstablyand 81 (a) (b) Figure4.2:ThetorsionalAlfvenwavepulse.Showninthesepanelsareplotsattime t =0 : 156 of(a)theenergyand(b)thethermalpressure.Thesolutionwasobtainedonameshwith N =800. thenumericalresultsat t =0 : 156areshowninFigures4.2and4.3with N =800.Shown intheareplotsoftheenergy,thethermalpressure, u y , u z , B y ,and B z .Itis observedthatourmethodsuccessfullycapturesthetwodiscontinuitiesandtheresultsare comparablewiththosein[9,49].However,smallbumpscanstillbeobservedaroundone ofthediscontinuitiesofboth u y and u z .Aspointedoutin[9],thisisbecausetheMHD solverintroducedtoomuchnumericaldissipationtokeepthepressurepositive.Theprimary reasonistheRiemannsolveraroundthediscontinuitiesisnotselectiveenough. 4.4.2Testcasesinmulti-D Inthissubsection,weconsiderseveral2Dand3Dexamplestodemonstratetheaccuracy andofourpositivity-preservingmulti-DMHDsolverintheCTframework.Inthe followingtests,weimplementfourth-orderWENO-CT2DandWENO-CT3Dschemesasthe MHDsolver,towhichweapplyourpositivity-preservinglimiter.Unlessotherwisestated, weuse Option2 forthemulti-Dsimulation. 82 (a) (b) (c) (d) Figure4.3:ThetorsionalAlfvenwavepulse.Showninthesepanelsareplotsattime t =0 : 156 of(a) u y ,(b) u z ,(c) B y and(d) B z .Thesolutionwasobtainedonameshwith N =800. 4.4.2.1SmoothvortextestinMHD Weconsiderthesmoothvortexproblemwithnon-zeromagnetictodemonstratethe schemecanmaintainthedesignedaccuracywithintheCTframework.Weconsideramod- ofthesmoothvortexproblemconsideredin[3,48,80].Theinitialconditionisa meanw ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )=(1 ; 1 ; 1 ; 0 ; 1 ; 0 ; 0 ; 0) ; (4.37) 83 withperturbationson u x , u y , B x , B y ,and p : ( u x ;u y )= 2 ˇ e 0 : 5(1 r 2 ) ( y;x ) ; ( B x ;B y )= 2 ˇ e 0 : 5(1 r 2 ) ( y;x ) ; p = y (1 r 2 ) 2 8 ˇ 2 e 1 r 2 : Themagneticpotentialisinitializedas A z = 2 ˇ e 0 : 5(1 r 2 ) : Here r 2 = x 2 + y 2 . Wesetthevortexstrength =5 : 389489439and = p 2 suchthatthelowestpressure inthecenterofthevortexis5 : 3 10 12 .Similarto[48],weusecomputationaldomain ( x;y ) 2 [ 10 ; 10] [ 10 ; 10]suchthattheerrorfromtheboundaryconditionswillnot theoverallconvergencestudy.Theperiodicboundaryconditionareusedonall sides.Becausefourth-orderCTstepsareused,theoverallschemeisfourth-orderaccuracy. The L 1 -errorsand L 1 -errorsofthevelocityandmagneticfor t =0 : 05areshownin Table4.1,inwhichonecanconcludetheproposedpositivity-preservingschemecanmaintain fourth-orderaccuracyasexpected.Weremarkthatnegativepressureisobservedonmeshes coarserthan320 320whenthelimiterisnotapplied. 4.4.2.2Smoothvortextestinhydrodynamics Inthepreviousexample,wedemonstratetheproposedschemecanattainthedesignedorder ofaccuracywhenthereisverylowpressureinthesolution.Inthisexample,wedemonstrate theproposedschemecanalsoattainthedesignedorderofthebaseMHDsolverbythe 84 Table4.1:Accuracytestofthe2DvortexevolutioninMHD.Shownarethe L 1 -errorsand L 1 -errorsattime t =0 : 05ofthedensityascomputedbythepositivity-preservingWENO- CT2Dschemeatvariousgridresolutions. Mesh u x u y L 1 -error order L 1 -error order L 1 -error order L 1 -error order 40 40 7.38e-04 - 1.79e-02 - 8.03e-04 - 1.94e-02 - 80 80 7.20e-05 3.35 4.33e-03 2.05 7.36e-05 3.45 5.22e-03 1.90 160 160 3.46e-06 4.38 1.92e-04 4.49 3.72e-06 4.31 2.18e-04 4.58 320 320 1.80e-07 4.27 1.42e-05 3.76 1.96e-07 4.25 1.64e-05 3.73 Mesh B x B y L 1 -error order L 1 -error order L 1 -error order L 1 -error order 40 40 1.02e-03 - 1.49e-02 - 1.04e-03 - 1.57e-02 - 80 80 7.73e-05 3.73 1.27e-03 3.56 7.73e-05 3.75 1.16e-03 3.77 160 160 4.75e-06 4.03 8.25e-05 3.94 4.74e-06 4.03 7.16e-05 4.01 320 320 2.85e-07 4.06 7.66e-06 3.43 2.84e-07 4.06 6.36e-06 3.49 2Dvortexevolutionproblemfromhydrodynamics,wherethereareverylowdensityand pressure.Theinitialconditionconsistsofameanw ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )=(1 ; 1 ; 1 ; 0 ; 1 ; 0 ; 0 ; 0) ; (4.38) andperturbationson u x , u y andthetemperature T = p=ˆ : ( u x ;u y )= 2 ˇ e 0 : 5(1 r 2 ) ( y;x ) ; T = ( 1) 2 8 ˇ 2 e 1 r 2 ; withnoperturbationintheentropy S = p=ˆ .Here r 2 = x 2 + y 2 .Inthiscase,welet =1 : 4. Thecomputationaldomainis( x;y ) 2 [ 5 ; 5] [ 5 ; 5]withperiodicboundaryconditionon allsides.Wesetthevortexstrength =10 : 0828suchthatthelowestdensityandlowest pressureinthecenterofthevortexare7 : 8 10 15 and1 : 7 10 20 respectively.Theexact 85 Table4.2:Accuracytestofthe2Dvortexevolutioninhydrodynamics.Shownarethe L 1 -errorsand L 1 -errorsattime t =0 : 05ofthedensityascomputedbythepositivity- preservingWENO-CT2DschemewithCTstepsturnedoThesolutionsconvergeat orderaccuracy. Mesh L 1 -error order L 1 -error order 40 40 6.22e-03 - 2.30e-01 - 80 80 5.64e-04 3.46 3.98e-02 2.53 160 160 1.09e-05 5.69 1.16e-03 5.10 320 320 1.79e-07 5.93 2.75e-05 5.40 640 640 4.92e-09 5.19 8.71e-07 4.98 solutiontothisproblemisjusttheconvectionofthevortexwiththemeanvelocity. Weremarkthatthemagneticisinitializedaszeroandtheexactmagneticwill bezeroforallthefuturetime.AlthoughthenumericalsolutionsbytheMHDsolverhas anonzeromagnetictheofthecomputedmagnetictothewholesystemis tcomparedtotheotherquantities.Asaresult,theofCTstepsto theresultsisnegligibleforthiscase.Wesolvethisproblemwiththepositivity-preserving WENO-CT2DschemewithCTstepsturnedThe L 1 -errorsand L 1 -errorsofthedensity at t =0 : 05areshowninTable4.2.Weclearlyobserveaconvergenceofthebase schemes.Weremarkthatnegativedensityorpressureisobservedduringthecomputation ofallthemeshesinthetablewhentheproposedlimiterisnotapplied. 4.4.2.3Rotatedvacuumshocktubeproblem Inthisexample,weconsiderthevacuumshocktubeproblemrotatedbyanangleof ina 2Ddomain.Theinitialconditionsinthiscaseare ( ˆ;u ? ;u k ;u z ;p;B ? ;B k ;B z )= 8 > > > < > > > : (10 12 ; 0 ; 0 ; 0 ; 10 12 ; 0 ; 0 ; 0)if ˘< 0 ; (1 ; 0 ; 0 ; 0 ; 0 : 5 ; 0 ; 1 ; 0)if ˘> 0 : (4.39) 86 (a) (b) (c) (d) Figure4.4:Rotatedvacuumshocktubeproblem.Showninthesepanelsareplotsattime t =0 : 1of(a) ˆ ,(b) ˆ cutaty=0,(c)thethermalpressureand(d)thethermalpressurecut aty=0.40equallyspacedcontoursareusedforthecontourplots.Thesolidlinesin(b) and(d)are1Dhighlyresolvedsolutions.Thesolutionwasobtainedona240 100mesh. where ˘ = x cos + y sin and = x sin + y cos ; (4.40) where u ? and B ? areperpendiculartotheshockinterface,and u k and B k areparallelto theshockinterface.Themagneticpotentialisinitializedas A z (0 ;˘ )= 8 > > > < > > > : 0if ˘ 0 ; ˘ if ˘ 0 : (4.41) 87 Wesolvethisproblembythepositivity-preservingWENO-CT2Dschemeonthecomputa- tionaldomain( x;y ) 2 [ 0 : 6 ; 0 : 6] [ 0 : 25 ; 0 : 25]witha240 100mesh. =tan 1 (0 : 5). Zero-orderextrapolationboundaryconditionsareusedontheleftandrightboundaries.On thetopandbottomboundaries,allthequantitiesaresettodescribetheexactmotionofthe shock. ThesolutionsareplottedinFigure4.4,wherethe1Dcutofdensityandpressureat y =0 isalsoplottedtocomparewiththe1Dhighlyresolvedresults.Weclearlyobservethatthe 2Dsolutionisconsistentwiththe1Dsolution.Withoutthelimiter,negativedensityand pressureareobservedinnumericalsolutions,whichquicklyleadstoblow-upofthenumerical simulation. 4.4.2.42Dblastproblem Intheblastwaveproblem,astrongfastmagnetosonicshockformulatesandpropagates intothelow- plasmabackground,whichwilllikelyleadtonegativedensityorpressure innumericalsolutions.Inthissubsection,werstinvestigatea2Dversionoftheproblem [8,10,49].Thecomputationaldomainis( x;y ) 2 [ 0 : 5 ; 0 : 5] [ 0 : 5 ; 0 : 5]withw boundaryconditionsonallthefoursides. Theinitialconditionsoftheproblemconsistofaninitialbackground: ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )=(1 ; 0 ; 0 ; 0 ; 0 : 1 ; 100 = p 8 ˇ; 100 = p 8 ˇ; 0) ; (4.42) andacircularpressurepulse p =1000withinaradius r =0 : 1fromthecenterofthedomain. 88 (a) (b) (c) (d) Figure4.5:2Dblastproblem.Showninthesepanelsareplotsattime t =0 : 01of(a) ˆ ,(b)thethermalpressure,(c)thenormof u and(d)themagneticpressure.40equally spacedcontoursareusedforeachplot.Thesolutionwasobtainedona256 256meshby positivity-preservingWENO-CTscheme. Theinitialscalarmagneticpotentialissimplygivenby A z =100 = p 8 ˇy 100 = p 8 ˇx: (4.43) Thesolutioniscomputedona256 256mesh.ShowninFigure4.5areplotsofthesolution. Thesolutionshowsgoodagreementwiththosein[8,49]. InTable4.3,weusethisexampletocomparefourtschemes,WENO-HCL, 89 (a) (b) (c) (d) Figure4.6:2Dblastproblem.Showninthesepanelsareplotsattime t =0 : 01of(a) ˆ ,(b) thethermalpressure,(c)thenormof u and(d)themagneticpressure.40equallyspaced contoursareusedforeachplot.Thesolutionwasobtainedona256 256meshbythe WENO-HCLscheme. WENO-CT-OP1,WENO-CT-OP2andPP-WENO-CT-OP2.HereWENO-HCLisreferred tothebaseWENO-HCLschemewithoutCTorthelimiter.WENO-CT-OP1isreferred totheWENO-CT2Dschemechoosing Option1 withoutthelimiter.WENO-CT-OP2is referredtotheWENO-CT2Dschemechoosing Option2 withoutthelimiter.Finally,PP- WENO-CT-OP2isreferredtothepositivity-preservingWENO-CT2Dschemewith Option 2 chosen.FromTable4.3,weobservethatthebaseWENO-HCLschemeisunstablefor 90 Table4.3:Comparisonsoftschemessolvingthe2Dblastproblem.Thecolumnof P indicateswhetherthenumericalsolutionsremainpositiveinthesimulations.Thecolumn of S indicateswhetherthesimulationsrunstablyto t =0 : 01.Inordertomakeafair comparison,thepositivityofdensityandpressureisonlycheckedateachtimestep t n . Mesh WENO-HCL WENO-CT-OP1 WENO-CT-OP2 PP-WENO-CT-OP2 P S P S P S P S 150 150 No No No No No Yes Yes Yes 200 200 No Yes No No No Yes Yes Yes 256 256 No Yes No No No No Yes Yes theresolution150 150andbecomesstableinthehigherresolutions.WENO-CT-OP1is unstableforeachresolutionandapplyingthepositivity-preservinglimiterwillnotbeableto stabilizethisbecausethenegativepressureisfromthecorrectionstepofthemagnetic WENO-CT-OP2isstableinthelowerresolutionbutbecomesunstablefortheresolution 256 256,andnegativepressureisobservedinalltheresolutions.Finally,thepositivity- preservingWENO-CTschemeisstableforalltheresolutions.Fromthoseresults,itisvery clearthatthepositivity-preservingWENO-CTschemeisthemoststablemethodinthese fourmethods. AnotherconcernforthisCTframeworkisthattheenergyisnotconservinginour positivity-preservingWENO-CTschemedueto Option2 .Wealsousethisexampleto studythisissue.WecomparetheresultsbythebaseWENO-HCLschemeandthepositivity- preservingWENO-CTscheme.InFigure4.6,weshowtheresultsbythebaseWENO- HCLschemewiththesameresolutionasFigure4.5.Theresultslooksimilartothoseby thepositivity-preservingWENO-CTscheme,excepttherearesomeunphysicaloscillations aroundthecenterregioninFigure4.6.Thatisduetothedivergenceerrorinthebasescheme. Ifweplotthedivergenceerrorinthetimedomain,wecanclearlyseethatthedivergence errorofthepositivity-preservingWENO-CTschemestaysaround10 12 ,whiletheerrorof 91 (a) (b) Figure4.7:ComparisonsbetweenWENO-HCLschemeandpositivity-preservingWENO-CT schemeforthe2Dblastproblem.Showninthesepanelsareplotsof(a)thedivergenceerror and(b)therelativeerrorofthetotalenergyinthewholedomain.The x -axisdenotesthe time t 2 [0 ; 0 : 01]. theWENO-HCLschemeisaround10 0 duringthesimulation.Herethedivergenceerroris as L 1 -normof r B ,wherethenumerical operatorisasaregularfourth- ordercentraldiscretization.AsacommondrawbackintheCTframework when Option2 ischosen,thecorrectionstepleadstoalossoftheconservationofthetotal energy. Tostudythisissue,wecomputethetotalenergiesinthewholedomainsolvedbytwo schemesandcomparethemwiththeinitialnumericalvalue.Theresultsareplottedin Figure4.7.Wecanseetherelativeerrorofthepositivity-preservingschemeisaround 10 3 ,whiletheconservativeWENO-HCLschemehasanerrorabout10 12 .Ontheother hand,whenweperformaconvergencestudyforthetotalenergyinthewholedomain, aconvergencehasbeenobservedforthepositivity-preservingschemewhenthe computedtotalenergiesarecomparedwiththeexactvalue,althoughtheschemeisnot energy-conserving.However,weremarkthattheconservationofenergyisimportantforsome problems,suchasthoseinvolvingnonlinearstrongdiscontinuities.Ahigh-orderpositivity- 92 (a) (b) Figure4.8:3Dblastproblem.Showninthesepanelsare3Dplotsattime t =0 : 01of(a) ˆ and(b)thethermalpressure.Thesolutionwasobtainedona150 150 150mesh. preservingconservativeschemewiththedivergenceerrorbeingcontrolledwillbepartofour futurework.However,itisverycult,ifnotimpossible,intheCTframeworktosatisfy alltherequirementssimultaneously.Abetterwaytocontrolthedivergenceerrorisneeded forthispurpose. 4.4.2.53Dblastproblem Thelastproblemweinvestigateisafully3Dversionoftheblastproblem.Itisusedto testthebehaviorofthepositivity-preservingWENO-CT3Dscheme.Theinitialconditions consistofaninitialbackground, ( ˆ;u x ;u y ;u z ;p;B x ;B y ;B z )=(1 ; 0 ; 0 ; 0 ; 0 : 1 ; 100 = p 8 ˇ; 100 = p 8 ˇ; 0) ; (4.44) andasphericalpressurepulse p =1000withinaradius r =0 : 1fromthecenterofthe domain.Theinitialconditionsforthemagneticpotentialare A (0 ;x;y;z )=(0 ; 0 ; 100 = p 8 ˇy 100 = p 8 ˇx ) : (4.45) 93 (a) (b) (c) (d) Figure4.9:3Dblastproblem.Showninthesepanelsareplotsattime t =0 : 01andcutat z =0of(a) ˆ ,(b)thethermalpressure,(c)thenormof u and(d)themagneticpressure.40 equallyspacedcontoursareusedforeachplot.Thesolutionwasobtainedona150 150 150 mesh. Thecomputationaldomainis[ 0 : 5 ; 0 : 5] 3 .wboundaryconditionsareusedonall sides.Thenumericalsimulationisperformedona150 150 150mesh.Todistinguish this3Dcasefromthe2Dblastcase,wepresentthe3Dplotsofthedensityandpressurein Figure4.8,whichclearlyindicatesitssphericalstructures.InFigure4.9weshowtheresults ofthesolutionscutat z =0.Thesolutioniscomparabletothe3Dresultsin[31,53,85]. Wenotethatnegativepressureisobservedattime t =0 : 0033ifthepositivity-preserving 94 limiterisnotapplied. 95 Chapter5 Conclusionandfuturework 5.1Conclusion InthisworkwedevelopedaclassofmethodsforsolvingtheidealMHD equations.Asummaryofthekeyfeaturesoftheproposednumericalmethodarelisted below: 1. Allquantities,includingallcomponentsofthemagneticandmagneticpotential, aretreatedaspointvaluesonthesamemesh(i.e.,thereisnomeshstaggering). 2. ThebaseschemeistheorderFD-WENOschemeofJiangandShu[39].Withthis methodweareabletoachievehigh-orderusingdimension-by-dimension operators,insteadofthemorecomplicatedspatialintegrationandmultidimensional reconstructionsusedbyHelzel etal. [35]. 3. Thecorrectedmagneticiscomputedviafourth-orderaccuratecentraldif- ferenceoperatorsthatapproximatethecurlofthemagneticvectorpotential.These operatorsarechosentoproduceacorrectedmagneticldthatexactlya discretedivergence-freecondition. 4. Alltime-steppingisdonewiththeSSP-RKmethods.InChapter3,weuseafourth- orderversionandinChapter4weuseathird-orderversion. 96 5. Usingaparticulargaugecondition,themagneticvectorpotentialismadetosatisfya weaklyhyperbolic,non-conservative,hyperbolicsystem.Thissystemissolvedusinga moversionoftheFD-WENOschemedevelopedforHamilton{Jacobiequations [38].For3Dproblems,speciallimitersbasedonialresistivityareintroducedto helpcontrolunphysicaloscillationsinthemagnetic 6. Apositivity-preservinglimiterisproposedandappliedbymodifyinghigh-orderWENO- HCLwiththerst-orderLax{Friedrichssoastoproducepositivedensityand pressure. InChapter3thenumericalmethodsweretestedonseveral2D,2.5D,and3Dtestprob- lems,allofwhichdemonstratetherobustnessofourapproach.Onsmoothproblems,we achievefourth-orderaccuracyinallcomponents,includingthemagneticandthemag- neticpotential.Forproblemswithshocks,weareabletoaccuratelycapturetheshockwaves withoutintroducingunphysicaloscillationsinanyofthesolutioncomponents.Inaddition, thecloud-shockinteractionproblemsalsoindicatedthatthereisapossibleadvantageof usingahigh-ordermethodcomparedtotraditionalsecond-ordermethods.Forinstance, usinga128 128meshinourmethods,weareabletoseethesamestructuresthatcanonly beobservedbyasecond-ordervolumemethodsonmuchgridresolutions.This phenomenonisobservedinboth2Dand3D.Anotheradvantageoftheproposedmethodsis thattheydonotinvolveanymultidimensionalreconstructionsinanystep,whiletraditional high-ordervolumemethodscommonlyneedseveralmultidimensionalreconstructions ineachgridcell.Forinstance,forthesameresolutionona3Dsimulation,our codeuseslessCPUtimethanthethird-ordervolumecodein[35]. 97 InChapter4wedemonstratedtheenessandoftheresultingpositivity- preservingschemeswith1D,2Dand3Dproblems.Throughsmoothproblems,weshowed theproposedlimiterwillmaintainthedesignedorderofaccuracyofthebaseschemes.We alsotestedtheschemeswithseveralpracticalproblemssuchasthelow- problemsandfound theproposedlimitercanincreasethenumericalstabilityoftheschemes.Inparticular,we alsostudiedtheofenergycorrectionstepsintheconstrainedtransportmethodbya blastwaveproblem. 5.2Futurework Thenumericalschemesasdevelopedsofarcanonlybeusedtosolveproblemsoneithera uniformgridoronasmoothlyvaryingmappedgrid,whichisacommondisadvantageofFD ENO/WENOschemes.Thus,ourmethodsarelesscomparedtothevolume CTmethodsdevelopedin[35],inwhichtheschemehasbeensuccessfullyextendedtonon- smoothlyvaryinggrids.However,apromisingapproachforovercomingthisrestrictionisto placetheexistingWENO-CTmethodforidealMHDintoanAMRframework.Sinceour methodsarefullyexplicitandfullyunstaggered,itispossibletoincorporatethemintothe WENO-AMRframeworkdevelopedbyShen etal. [66].Tobetterachievethegoal,currently weareworkingonahigh-ordersingle-stepsingle-stageidealMHDsolverbyextendingthe workin[63]. InadditiontoAMRapplication,wearealsointerestedinincludingnon-idealterms,such asHalltermandresistivityterm,intotheidealMHDsolversothatitcansimulatemore physicalphenomenasuchasmagneticreconnection.Thenon-idealtermwillpotentially becomeveryifitsistforthesystem.Consequently,thistermneeds 98 tobetreatedimplicitlyandonemayneedothertimesteppingapproachessuchasimplicit- explicittimeintegrationmethodsandfullyimplicitmethods.Recently,wehavebeenalso applyingtheproposedschemestoseveralpracticalproblems,suchassimulating Z-pinchplasmas. 99 APPENDIX 100 WENOreconstruction ThemainideaofWENOreconstructionistocomputeastencilusinga weightedaverageofseveralsmallerstencils.Theweightsarechosenbasedonthesmoothness ofthesolutiononeachofthesmallerstencils.Thefullprocedurecanbefoundin[39, 64,67,68].Forcompleteness,wealsoincludeabriefdescriptionofthe5 th -orderWENO reconstructionasusedinthisworkandtheoperator WENO5 thatwasusedinSection 3.1. Weconsidertheproblemonauniformgridwith N +1gridpoints, a = x 1 2