HIGHORDERFINITEDIFFERENCEWENOSCHEMESFORIDEALMAGNETOHYDRODYNAMICSByXiaoFengADISSERTATIONSubmittedtoMichiganStateUniversityinpartialentoftherequirementsforthedegreeofMathematics|DoctorofPhilosophy2017ABSTRACTHIGHORDERFINITEDIFFERENCEWENOSCHEMESFORIDEALMAGNETOHYDRODYNAMICSByXiaoFengInthisdissertationweproposetwohighordernumericalschemesforsolvingtheidealmagnetohydrodynamic(MHD)equations.Theschemeissingle-stagesingle-step,maintainsadivergence-freeconditiononthemagneticandhasthecapacitytopreservethepositivityofthedensityandpressure.Toaccomplishthis,weuseaTaylordiscretizationofthePicardintegralformulation(PIF)oftheWENOmethodproposedin[SINUM,53(2015),pp.1833{1856].Weusetheversionwherexesareexpandedtothird-orderaccuracyintime,andforthevariablesspaceisdiscretizedusingtheclassicaleWENOdiscretization.Weuseconstrainedtransportinordertoobtaindivergence-freemagneticwhichmeansthatwesimultaneouslyevolvethemagnetohydrodynamic(thathasanevolutionequationforthemagneticandmagneticpotentialequationsalongsideeachother,andsetthemagnetictobethe(discrete)curlofthemagneticpotentialaftereachtimestep.Inthiswork,wecomputethesederivativestofourth-orderaccuracy.Inordertoretainasingle-stage,single-stepmethod,wedevelopanovelLax-Wdiscretizationfortheevolutionofthemagneticpotential,wherewestartwithtechnologyusedforHamilton-Jacobiequationsinordertoconstructanon-oscillatorymagneticTheendresultisanalgorithmthatissimilartoourpreviouswork[JCP,268,(2014),pp.302{325],butthistimethetimesteppingisreplacedthroughaTaylormethodwiththeadditionofapositivity-preservinglimiter.Finally,positivitypreservationisrealizedbyintroducingaparameterizedlimiterthatconsidersalinearcombinationofhighandlow-ordernumericalThechoiceofthefreeparameteristhengiveninsuchawaythatthearelimitedtowardsthelow-ordersolveruntilpositivityisattained.Giventhelackofadditionaldegreesoffreedominthesystem,thispositivitylimiterlacksenergyconservationwherethelimiterturnson.ThesecondschemeisbasedonanalternativeformulationoftheWENOscheme[JCP,77(1988),pp.439{471,SIAMJ.Sci.Comput.,35(2013),pp.A1137{A1160,andMethodsAppl.Anal.,21(2014),pp.1{30],whichcanbeappliedoncurvilinearmeshes.ThisschemecomputesthenumericalbyaTaylorexpansioninspace.ThelowestordertermintheexpansioniscomputedbyapplyingaRiemannsolvertoone-sidedapproximationstoconservedquantitiesobtainedfromWENOinterpolation.Thehigherordertermsintheexpansionarecomputedbycentrales.AnadditionallimiterbasedonthesmoothnessindicatorsfromtheWENOinterpolationisappliedtocontroltheoscillationsinthehigherorderterms.Constrainedtransportisdoneinawaysimilarto[JCP,268,(2014),pp.302{325].Apositivitypreservinglimitersimilartotheoneintheschemeisappliedinordertoenhancetherobustnessofthescheme.WhentheRiemannsolverusedinthisschemeischosentobeahighresolutionsolversuchastheHLLDsolver,theschemeexhibitsbetterresolutionthanschemesbasedonWENOreconstructionofandLax-Friedrichssplitting.WepresentnumericalresultsforseveralstandardtestproblemsincludingsmoothAlfvenwaveproblems(toverifyformalorderofaccuracy),shocktubeproblems(totesttheshock-capturingabilityofthescheme),Orszag-Tang,andcloudshockinteractions.Theseresultsasserttherobustnessandverifythehigh-orderofaccuracyoftheproposedscheme.ACKNOWLEDGMENTSIwouldliketoexpressmysinceregratitudetomyadvisor,ProfessorAndrewChristlieb,forhiscontinuoussupportofmyPhDstudy.Itwastrulyjoyfulandrewardingtoworkunderhisguidance.Iwouldalsoliketothanktherestofmydissertationcommittee:Profes-sorYingdaCheng,ProfessorSeanCouch,ProfessorKeithPromislow,andProfessorJianliangQian,fortheirinsightfulcommentsandadvice.Mysincerethanksalsogoestomyfellowgroupmates,friends,andcollaborators,Dr.DavidSeal,Dr.QiTang,andDr.YanJiang,fornumerousstimulatingdiscussions.IwouldalsoliketothankMs.BarbaraMillerfortheexampleofprofessionalismshesetforme,andforhertremendoushelp.Iamextremelygratefulformywife,ZhijieYu,forsheisthebestlifecompanionIcanpossiblyimagine.Thankyouandlet'scontinueourexplorationstogether!Lastbutnottheleast,Iwouldliketothankmyparentsfortheirsupportthroughouttheyears.Thankyou!ivTABLEOFCONTENTSLISTOFTABLES....................................viiLISTOFFIGURES...................................viiiChapter1Introduction...............................1Chapter2TheidealMHDequations.......................92.1HyperbolicityoftheidealMHDequations...................102.2DiscontinuitiesintheidealMHDequations...................12Chapter3TheWENOmethodology.......................163.1WENOinterpolation...............................163.2WENOreconstruction..............................20Chapter4Ahigh-orderpositivity-preservingsingle-stagesingle-stepmethodfortheidealmagnetohydrodynamicequation...........224.1Aconstrainedtransportframework.......................224.2TheTaylor-discretizationPIF-WENOmethod.................264.3Theevolutionofthemagneticpotentialequation...............274.3.1The2Dmagneticpotentialequation..................294.3.2The3Dmagneticpotentialequation..................354.4Positivitypreservation..............................384.5Numericalresults.................................454.5.1SmoothAlfvenwave...........................454.5.1.1SmoothAlfvenwave:The2Dproblem............464.5.1.2SmoothAlfvenwave:The3Dproblem............464.5.22Drotatedshocktubeproblem.....................484.5.32DOrszag-Tangvortex..........................534.5.42Drotorproblem.............................554.5.5Cloud-shockinteraction.........................594.5.5.1Cloud-shockinteraction:The2Dproblem..........594.5.5.2Cloud-shockinteraction:The3Dproblem..........614.5.6Blastwaveexample............................634.5.6.1Blastwaveexample:The2Dproblem............634.5.6.2Blastwaveexample:The3Dproblem............644.5.7Errorsinenergyconservation......................64Chapter5Ahigh-orderschemeforidealmagnetohydro-dynamicsoncurvilinearmeshesbasedonanalternativefor-mulationoftheWENOschemeandconstrainedtransport..705.1AnalternativeformulationoftheWENOscheme.............71v5.1.1Introduction................................715.1.2Alimiteronthehigherordertermsinthenumerical.......735.1.3One-dimensionalsystems.........................755.1.4Multidimensionalcases..........................785.1.5Curvilinearmeshes............................795.2HLLtypeRiemannsolversforMHDequations.................805.2.1TheHLLapproximateRiemannsolver.................825.2.2TheHLLCapproximateRiemannsolver................835.2.3TheHLLDapproximateRiemannsolver................875.3Constrainedtransport..............................945.4Apositivity-preservinglimiter..........................965.5Numericalresults.................................975.5.12DsmoothAlfvenwaveproblem....................975.5.2Brio-Wushocktube...........................985.5.2.11Dshocktube.........................995.5.2.22Drotatedshocktube.....................1025.5.32DOrszag-Tangvortex..........................1025.5.42Dcloud-shockinteraction........................108Chapter6Conclusionsandfuturework.....................113BIBLIOGRAPHY....................................116viLISTOFTABLESTable4.1:2DsmoothAlfvenwave.Here,weshowL1-errorsatatimeoft=1:0forthesolutionwithandwithouttheenergy\correction"step.Table(a)hastheenergycorrectionturnedandTable(b)hastheenergycorrectionturnedon.Becausetimeisonlydiscretizedtothird-orderaccuracy,weobservethepredictedthird-orderaccuracyofthesolverhere........47Table4.2:2DsmoothAlfvenwave.Here,weshowL1-errorsatashorttimeoft=0:01forthesolutionwithandwithouttheenergy\correction"step.Table(a)hastheenergycorrectionturnedandTable(b)hastheenergycorrectionturnedon.Here,weetfasterthanxinordertoexposethespatialorderofaccuracyofthesolver.Becauseweonlyuseafourth-orderaccuratespatialdiscretizationforrA,weonlyobservefourth-orderaccuracydespitethefactthatthevariablesarediscretizedtoorderaccuracy.................................47Table4.3:3DsmoothAlfvenwave.InthistableweshowtheL1-errorsatamoderatetimeoft=1:0.InTable(a)thepositivity-preservinglimiterisandthethepositivity-preservinglimiter(andtheenergycorrectionstep)isturnedonfortheresultsinTable(b).Becausetimeisdiscretizedtothird-orderaccuracy,themethodisformallyonlythird-orderaccurateintime.InTable4.4werunthesolvertoashorttimeinordertoexposethespatialorderofaccuracy............................48Table4.4:3DsmoothAlfvenwave.Here,weshowtheL1-errorsatashorttimeoft=0:01.Inaddition,wetfasterthanthemeshspacinginordertoextractthespatialorderofaccuracy.Table(a)hasthepositivity-preservinglimiterturnedandTable(b)hasthepositivity-preservinglimiter(aswellasthecorrectionstep)turnedon.Theresultsarealmostidentical.....................................49Table5.1:L1-errorsin2DsmoothAlfvenproblem.AlternativeWENOwithLax-FriedrichsConstrainedtransportandpositivity-preservinglimiteron.98Table5.2:L1-errorsin2DsmoothAlfvenproblem.AlternativeWENOwithHLLDConstrainedtransportandpositivity-preservinglimiteron......99viiLISTOFFIGURESFigure4.1:Stencilneededformagneticpotential.ThevaluesofAzateach\"(savethecentermostvalue)isneededtocorrectthevalueofBatthecenter.TheshadedpointsareneededinEqns.(4.26)and(4.30){(4.34)tocomputethevalueofAzattherightmost\"point.Notethatthisstenciliscontainedinthestencilshownin[20],Fig.2...........34Figure4.2:Therotatedshocktubeproblem.Here,wecomparethesolverwith(leftpanels)andwithout(rightpanels)constrainedtransportturnedon.Auniformmeshofsize200100isusedforbothsimulations.Theangleofrotationfortheinitialconditionsis˚=tan1(0:5),and30equallyspacedcontoursrangingfromtheminumumtothemaximumofeachfunctionareusedforthetoptwopanels.ThecontourplotforthesolutionwithoutCTcontainssmallwigglesthataremuchmorepronouncedwhenslicesofthesolutionaresampled.Tothisend,asliceofthesolutionalongy=0ispresentedinthebottomtwopanels.FurtherevolutionproducesasolutionthatcausesthecodetofailinthecasewhereCTisturnedThepositivity-preservinglimiteristurnedinordertoexercisethecode.Thesolidlinesinthebottomimagesarereferencesolutionsthatarecomputedbysolvingtheequivalent1Dshockproblemonauniformmeshwith50,000pointswiththeWENOmethod.51Figure4.3:2Drotatedshocktubeproblem.Theleftpanelshaveconstrainedtrans-portturnedon,andtherightpanelshaveconstrainedtransportturnedComponentsofthemagneticatt=0:2alongtheslicey=0.Themeshsizeis200100.Thepositivity-preservinglimiteristurnedinordertoexercisethecode.EachsolidlineisthereferencesolutiondescribedinFig.4.2.Weobservethatconstrainedtransport(withtheHamilton-Jacobisolver)allowsustonumericallycomputemagneticwithfarfeweroscillationsthanwouldotherwisebeobtainable......52Figure4.4:TheOrszag-Tangproblem.Weshowdensitycontourplotsat(a)t=0:5,(b)t=2,(c)t=3,(d)t=4.ThesolutioniscomputedusingPIF-WENOwithconstrainedtransportonandpositivity-preservinglimiteron,ona192192mesh.Atotalof15equallyspacedcontoursareusedforeachgraph......................................54viiiFigure4.5:TheOrszag-Tangproblem.Densityplotsatt=1:5forsolutionscomputedusingtinthenumericalscheme.(a)PIF-WENO,withconstrainedtransportturned(b)PIF-WENOwithconstrainedtransportandpositivitypreservinglimiterturnedon;(c)PIF-WENOwithconstrainedtransportturnedonandpositivitypreservinglimiterturned(d)Theslicealongy=0:5072.Thesolutionsarecomputedwitha192192mesh.Atotalof15equallyspacedcontoursareusedforeachof(a),(b),and(c)...............................56Figure4.6:TheOrszag-Tangproblem.Pressureplotsatt=3:0forsolutionscom-putedusingtinthenumericalscheme.(a)PIF-WENOwithconstrainedtransportturnedonandpositivitypreservinglimiterturned;(b)PIF-WENOwithconstrainedtransportandposi-tivitypreservinglimiterturnedon;(c)Theslicealongy=1:9799.Thesolutionsarecomputedwitha192192mesh.Atotalof15equallyspacedcontoursareusedforeachof(a)and(b)...............57Figure4.7:The2Drotorproblem.Shownhereare(a)thepseudocolorplotofthemagneticpressureand(b)themagneticline.Thissolutioniscom-putedonameshofsize400400.Constrainedtransportandpositivity-preservinglimitersareturnedon.......................58Figure4.8:The2Dcloud-shockinteractionproblem.Here,werunthesolvertoatimeoft=0:06.Inthethreepanels,weshowSchlierenplotsfor(a)thenaturallogofthedensity,(b)thenormofthemagneticand(c)thepressureforameshofsize256256.Thesameresultsforameshofsize512512arepresentedinpanels(d){(f),whereweobservemuchhigherresolutionfortheproblem.Constrainedtransportandpositivity-preservinglimiterareturnedonforbothsimulations............60Figure4.9:The3Dcloud-shockinteractionproblem.Schlierenplotsofln(ˆ).Thesolutionhereiscomputedusinga256256256mesh.Cross-sectionsaty=0:5andz=0:5fortheregiony0:5andz0:5areshown.Constrainedtransportandpositivity-preservinglimiterareturnedon...62Figure4.10:2Dblastproblem.Shownherearethecontourplotsatt=0:01of(a)density,(b)thermalpressure,(c)magnitudeofvelocity,and(d)magneticpressure.Atotalof40equallyspacedcontoursrangingfromthemintothemaxofthefunctionareusedforeachplot.Themeshsizeis256256.65Figure4.11:3Dblastproblem.Shownherearethepseudocolorplotsatt=0:01of(a)densityand(b)pressure.Themeshsizeis150150150.InFigure4.12weplotacutofthesolutionalongz=0.Thepositivity-preservinglimiterisrequiredtosimulatethisproblem.....................66ixFigure4.12:3Dblastproblem.Shownarethecontourplotsattimet=0:01cutatz=0of(a)density,(b)thermalpressure,(c)normofvelocity,and(d)magneticpressure.Thesolutionisobtainedusinga150150150mesh.Atotalof40equallyspacedcontoursareusedforeachplot........67Figure4.13:EnergyconservationerrorsforthesmoothAlfventestcases.Shownhereareresultsforthe(a)2DsmoothAlfven,and(b)3DsmoothAlfventestcases.Inordertoextracttheerrors,weplottheresultsonasemi-logscalebecauseotherwisetheresultsareindiscerniblefromthet-axis.Forthissmoothtestcase,theofthetheenergycorrectionstep(andhencethepositivity-preservinglimiter)isnegligiblebecausethesolutionremainssmoothfortheentiresimulation...................68Figure4.14:Energyconservationerrors.Shownhereareconservationerrorswhenthepositivity-preservinglimiter(andhencetheenergycorrectionstep)isturnedon.(a)2DOrszag-Tangproblem,(b)2Drotorproblem,(c)2Dblastproblem,and(d)3Dblastproblem.Notethattherotorandblastproblemsrequiretheapplicationofapositivity-preservinglimiterinordertorun,butthiscomesattheexpenseoflosingenergyconservation.FortheOrszag-Tangtestproblem,weruntoalatetime.Again,weobservethattheerrorsinenergyconservationdecreaseasthemeshis....................................69Figure4.15:Energyconservationerrors.Here,weshowresultsforconservationerrorswhenthepositivity-preservinglimiter(andhencetheenergycorrectionstep)isturnedThesolveranalyticallyconservesthediscretetotalenergyuptomachineprecision.Shownhereareresultsforthe(a)2DsmoothAlfven,and(b)3DsmoothAlfventestcases.Notethelogarith-micscalefortheaxes,andthattheseerrorsarenumericallyconserveduptomachineprecision..............................69Figure5.1:Aschematicdiagramofthephysicalandcomputationaldomainforcurvi-linearmeshes..................................79Figure5.2:Perturbedmeshusedfor2DAlfvenproblem.3232mesh.(˘;)2[0;1]2.x=˘+0:01sin(2ˇ2),y=+0:02sin(2ˇ˘4)..............98Figure5.3:Brio-Wushocktube.t=0:2.AlternativeWENOwithLax-FriedrichsAuniformmeshof200pointswereused.Positivity-preservinglim-iteristurnedon................................100Figure5.4:Brio-Wushocktube.t=0:2.AlternativeWENOwithHLLDAuniformmeshof200pointswereused.Positivity-preservinglimiteristurnedon....................................100xFigure5.5:Brio-Wushocktube.t=0:2.AlternativeWENOwithLax-FriedrichsAnon-uniformmeshof200pointswereused.Positivity-preservinglimiteristurnedon..............................101Figure5.6:Brio-Wushocktube.t=0:2.AlternativeWENOwithHLLDAnon-uniformmeshof200pointswereused.Positivity-preservinglimiteristurnedon..................................101Figure5.7:RotatedBrio-Wushocktube.t=0:2.AlternativeWENOwithLax-Friedrichs200100uniformmesh.Constrainedtransportandpositivity-preservinglimiterturned...................103Figure5.8:RotatedBrio-Wushocktube.t=0:2.AlternativeWENOwithLax-Friedrichs200100uniformmesh.Constrainedtransportandpositivity-preservinglimiterturnedon....................104Figure5.9:RotatedBrio-Wushocktube.t=0:2.AlternativeWENOwithHLLD200100uniformmesh.Constrainedtransportandpositivity-preservinglimiterturnedon..........................105Figure5.10:Orszag-Tangtestproblem.Contourplotofdensityatt=3.Aperturbed192192mesh.15equallyspacedcontourlines.AlternativeWENOwithLax-FriedrichsConstrainedtransportandpositivity-preservinglimiterturnedon................................106Figure5.11:Orszag-Tangtestproblem.Plotofdensityatt=3.Aperturbed192192mesh.15equallyspacedcontourlines.AlternativeWENOwithHLLDConstrainedtransportandpositivity-preservinglimiterturnedon..107Figure5.12:2Dcloud-shockinteraction.AlternativeWENOwithLax-Friedrichs256256uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon.................109Figure5.13:2Dcloud-shockinteraction.AlternativeWENOwithHLLD256256uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon....................109Figure5.14:Sketchofthenon-uniformmeshusedincloud-shockproblem.Actualmeshis256256.A3232meshisshownhereforclarity........110Figure5.15:2Dcloud-shockinteraction.AlternativeWENOwithLax-Friedrichs256256non-uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon.Thesecondrowshowstheregionnearthecenterofthecomputationaldomain.......111xiFigure5.16:2Dcloud-shockinteraction.AlternativeWENOwithHLLD256256non-uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon.Thesecondrowshowstheregionnearthecenterofthecomputationaldomain............112xiiChapter1IntroductionAplasmaisasubstancewhereatlylargeportionconsistsoffreechargedparticles,sothatduetothetheelectromagneticforces,whicharelong-rangeinnature,thesubstanceasawholeexhibitscollectivebehaviors[11].Magnetohydrodynamics(MHD)aremodelsforplasmaswhereitisassumedthattheplasmaisquasi-neutral(zeronetchargeinanyspatialscalesunderconsideration).TheidealMHDequationsfurtherassumeperfectconductivity(zeroelectricresistance).Althoughitisoneofthesimplestmodelsofplasmas,itisnonethelessa\remarkablyaccuratedescriptionofthelow-frequency,long-wavelengthdynamicsofrealplasmas"[72].Asaconsequence,theMHDequations,bothidealandnon-ideal,applicationsinawiderangeofareasincludingspacephysicsandcontrolledthermonuclearfusion[53,11].Inthisdissertation,weproposeandinvestigatetwohigh-ordernumer-icalschemesforidealMHDequationsbasedontheWeightedEssentiallyNon-Oscillatory(WENO)methodology.Theschemeisaasingle-stage,single-stepschemethatisde-signedwiththegoalofembeddinginAdaptiveMeshtframeworksandachievinglowstorageinmind.ThesecondschemeisaschemebasedonanalternativeformulationoftheWENOscheme,whichcanbeappliedoncurvilinearmeshes.Inbothschemes,weuseunstaggeredconstrainedtransport,whereamagneticpotentialisevolvedalongsidewiththeMHDquantities,tocontrolthedivergenceerrorofthemagnetic,andweusea1limitertoretainthepositivityofthedensityandpressureinMHD.High-ordernumericalschemesbasedontheEssentiallyNon-Oscillatory(ENO)[61]aswellastheWeightedENO(WENO)reconstructiontechnique[50,59,2,75,4,23,17]havebeenappliedsuccessfullytoidealMHDinthepasttwodecades.Thesehighorderschemesarecapableofresolvingcomplexfeaturessuchasshocksandturbulencesusingfewergridpointsthanlow-orderschemesforthesameleveloferror,asiscommonwithmanyhigh-ordershock-capturingschemes.Itoftenhappensinlarge-scaleMHDsimulationsthatthecomplexfeaturesareconcen-tratedinasmallportionofthesimulationdomain.Adaptivemesht(AMR)isatechniquethatisdesignedfortreatingsuchlocalityofcomplexityinhydrodynamicsandmagnetohydrodynamics.OneofthemainwithimplementinghighorderschemeswithinanAMRframeworkisthatboundaryconditionsfortheregionneedtobespinaconsistentmanner[62].ThisbecomesformultistageRunge-Kutta(RK)methods,becausehighordersolutionscannotbefoundifonesimplyuseshigh-orderinterpolatedvalues(intime)attheghostpointsthatarerequiredfortheintermediatestagesofthemethod.PreliminaryworkthatcombinesWENOspatialdiscretizationswithstrongstabilitypreservingRunge-Kutta(SSP-RK)time-steppingisconductedin[74,92],andveryrecentworkmakesuseofcurvilineargridstoextendmethodstoproblemswithgeometry[17].However,theauthorsin[17]useglobaltimesteps(whichprecludesthepossibilityofintroducinglocaltimestepping),andperhapsmoretroublesome,theydropmassconservationfortheirframeworktowork.InchoosingbuildingblocksforAMRcode,ithasbeenarguedthatsingle-stage,single-stepmethodsareadvantageous[24,5],partlybecausefewersynchronizationsareneededperstepthanmultistageRKmethods.Thefactthatsingle-stage,single-stepmethodsdonot2haveanissuewiththesesynchronizationsispossiblyoneofthereasonstheyhavegainedmuchattentioninthepasttwodecades,andonereasonwechoosetopursuethesemethods.Broadlyconstrued,thesemethodsarebasedonLaxandW'soriginalideaofusingtheCauchy-Kovalevskayaproceduretoconverttemporalderivativesintospatialderivativesinordertoanumericalmethod[54].Notablehigh-ordersingle-stage,single-stepmethodsincludetheArbitraryDERivative(ADER)methods[85,87],theLax-WWENOmethods[68],theLax-WdiscontinuousGalerkin(DG)meth-ods[67]andspace-timeschemesapplieddirectlytosecond-orderwaveequations[45,9].Ofthethreeclassesofhigh-ordermethodsbaseduponLax-Wtimestepping,onlytheADERmethodshavebeenappliedtomagnetohydrodynamics[5,4,13],whereassimi-larinvestigationshavenotbeendonefortheotherclasses.Anadditionaladvantagethatsingle-stagesingle-stepTaylormethodsistheirlow-storageopportunities.Thisrequirescare,becausethesemethodscaneasilyenduprequiringthesameamountofstorageastheirequivalentRKcounterpart(e.g.,ifeachtimederivativeisstoredinordertoreducecodingcomplexity).Thesingle-stagesingle-stepschemeweproposeinthisdissertationisbasedontheTaylordiscretizationofthePicardintegralformulationoftheWENO(PIF-WENO)method[20].ComparedwithotherWENOmethodsthatuseLax-Wtimediscretiza-tions[68],ourmethodhastheadvantagethatitsfocusisonconstructinghigh-orderTaylorexpansionsofthe(whichareusedtoaconservativemethodthroughWENOre-construction)asopposedtotheconservedvariables.Thisallows,forexample,theadaptationofapositivity-preservinglimiter,whichwedescribeinthisdocument.Theotherschemeweproposeinthisdissertationisbasedonanalternativefor-mulationofWENOschemes[51,52].Inthisalternativeformulation,thenumericalis3computedbyapplyingaRiemannsolvertoapproximationsofconservedquantitiesobtainedbyWENOinterpolation.ThisapproachcanbefoundintheearlyworkofShuandOsheronENOschemes[78].However,directreconstructionofnumericalhasenjoyedawidespreadpopularityfromearlyon(see[41]andthereferencestherein),inpartduetothelowercomputationalcost.Onedisadvantageofthisdirectreconstructionapproach,though,isthattheincorporationofupwindingintotheschemeiscommonlydonethroughtheveryeLax-Friedrichssplitting.ThealternativeformulationcanmitigatethisproblembyusinglesseRiemannsolvers.Formorediscussionontheadvantagesofthisalternativeformulation,wereferthereaderto[51,52].TheideaofusingRiemannsolversinnumericalschemesforhyperbolicconservationlawshasitsrootintheseminalworkofGodunov[35].Sincethen,atremendousamountofhasbeenputinpursuitofgoodapproximateRiemannsolvers.NotableworkonEulerequationsforhydrodynamicsincludetheworkofGodunov[36],Roe[69],Hartenetal[42],andToroetal[88].TheHarten-Lax-vanLeer(HLL)solverin[42]assumesoneintermediatestateintheapproximatesolutiontotheRiemannproblem.TheHLLC(Cforcontactdiscontinuity)solverin[88]improveduponthatbyassumingtwointermediatestatesandmakinguseofknowledgesptotheEulerequationtoachieveexactresolutionofisolatedcontactdiscontinuities.SimilarideaswereappliedtomagnetohydrodynamicsintheworkofGurski[40],Li[58],andMiyoshiandKusano[64].Ofthesesolvers,theHLLD(Dfordiscontinuities)solverin[64]assumesthelargestnumberofintermediatestates(four)intheapproximatesolutionandiscapableofresolvingcontact,rotational,andtangentialdiscontinuitiesexactly,and,whentheestimatesoffastestsignalspeedsareexact,alsoresolvingfastshocksexactly.ThissolverhassincebeenusedinmanymatureMHDsimulationcodesincludingFLASH[29]andAthena[82].4Weremarkherethattherehasbeenanongoing,inconclusivedebateonthecorrectconditionstosingleoutthephysicallycorrectweaksolutiontoaRiemannproblemforidealMHD.TheentropyconditionsthatworkedverywellforEulerequations(see,e.g.,[55])cannotservethispurposeforMHDequationsbecauseoflackofconvexityintheMHD[89].InnumericaltestsfortheshocktubeproblemexaminedbyBrioandWuin[15],intermediateshocks,whichweredeemedinadmissiblebyearliercriteria[48],showupubiquitouslyintheresults.Duetolackofdecisiveevidencefromphysicalobservationsandexperiments(however,see[32]),currentlythereisnoagreementonthephysicalrealityofintermediateshocks.TakahashiandYamada[83]givesanexcellentreviewofthecurrentcontroversyandalsoreportstheirownstudyofthisissue,whichareinfavorofFalleandKomissarov'sargument[31]thattheapparentintermediateshocksareactuallyaconsequenceofacertainsymmetryinthesetupoftheshocktubeproblem.Wealsoreferthereadertothesameauthors'work[84]forthedescriptionofanexactMHDRiemannsolver,whereonecouldswitchbetweenthevariousadmissibilityconditions.Anotherissuethatweaddressinoursecondschemeisthehandlingofcurvilinearmeshes(alsocalledstructuredgridsbysomeauthors).OnereasonforsolvingtheMHDequationsoncurvilinearmeshesisbecausethegeometryforcertainMHDapplicationsisbestdescribedbyanon-Cartesiancoordinatesystem[28,70,37,81].Morerecentworkinclude,forexample,thatofChacon[16]andStegmeiretal[80].Seealso[17]foraveryrecentpieceofworkonhighordernumericalschemesforgenerichyperbolicconservationlawoncurvilinearmeshes.Thewaywehandlecurvilinearmeshesinthe(second)schemeinthisdissertationfollowsthemethodin[52],whereahyperbolicsystemofconservationlawsontheCartesiancoordinatesistransformedtooneonthecurvilinearcoordinates.Wereferthereadertothereferencesthereinforfurtherdiscussions.5AnimportantissueinsimulationsofMHDsystemsisthecontrollingofthedivergenceerrorofthemagneticsincenumericalschemesbasedonthetransportequationsalonewill,ingeneral,accumulateerrorsinthedivergenceofthemagneticFailuretoaddressthisissuecreatesanunphysicalforceparalleltothemagnetic[14],andifthisisnottakencareof,itwilloftenleadtofailureofthesimulationcode.Populartechniquesusedtosolvethisprobleminclude(1)thenon-conservativeeight-wavemethod[38],(2)thepro-jectionmethod[14],(3)thehyperbolicdivergencecleaningmethod[27],and(4)thevariousconstrainedtransportmethods[1,8,23,25,30,33,43,44,71].othconductedanextensivesurveyin[90].BothschemesproposedinthisdissertationusestheunstaggeredconstrainedtransportframeworkproposedbyRossmanith[71].Thisframeworkevolvesavectorpotentialthatsitsonthesamemeshastheconservedquantities.ThisvectorpotentialisthenusedtocorrectthemagneticHistorically,theterm\constrainedtransport"hasbeenusedtorefertoaclassofmethodsthatincorporatethedivergence-freeconditionintothediscretizationofthetransportequationofthemagneticoftendoneinawaythatcanbeinterpretedasmaintaininganelectriconastaggeredmesh[30,90].Someauthorsactuallystilldistin-guishbetweenthistypeof\constrainedtransport"andthe\vectorpotential"approach[47].However,asiswell-known,evolvingavectorpotentialisconceptuallyequivalenttoevolv-inganelectricTheunstaggeredapproachhastheaddedbofeaseforpotentialembeddinginanAMRframework.Animportantpieceinanyvector-potentialbasedconstrainedtransportmethodisthediscretizationoftheevolutionequationofthevectorpotential.Thisevolutionequationisanonconservativeweaklyhyperbolicsystem,anditcanbetreatednumericallyfromthisviewpoint[44].AnalternativeapproachistoviewitasamodsystemofHamilton-Jacobi6equations[23].Bothschemesproposedinthisdissertationadoptsthelatterapproach.Thesingle-stage,single-stepschemeusesamethodinspiredbyaLax-WnumericalschemeforHamilton-Jacobiequationsthatwasproposedin[66].Theresistivitytermsusedin[23]areadaptedintothepresentsingle-stage,single-stepscheme.Thesecondscheme,whichisbasedonanalternativeformulationoftheWENOscheme,usesthesamemethodasin[23]toevolvethevectorpotential,withnecessarytransformationswhenappliedoncurvilinearmeshes.OnefurtherchallengefornumericalsimulationsofMHDisthatofretainingthepositivityofthedensityandpressure.Thisiscriticalwhenthethermalpressuretakesuponlyasmallportionofthetotalenergy,i.e.,whenthevalueissmall.Almostallpositivity-preservingmethodsexploitthepresumedpositivity-preservingpropertyofcertainlow-orderschemes,suchastheLax-Friedrichsscheme.Whereasearliermethodsoftenrelyonswitchingbetweenhigh-andlow-orderupdates[7,46],morerecentworktendstousecombinationsofthetwo.ExamplesincludetheworkofBalsara[3]andChengetal.[18],whichlimittheconservedquantitiesatcertainnodalpoints,andtheworkofChristliebetal.[22],whichcombineshigh-andlow-orderthroughasinglefreeparameterateachinterface.Inthecurrentwork,weadoptanapproachsimilartothatusedin[22].Thisapproachseeksasuitableconvexcombinationofthehigh-andlow-orderateachcellateachtimestep.Thislimiterweadoptisanadaptationofthemaximumprinciplepreserving(MPP)limiter[93](thathasitsrootsincorrectedtransportschemes[12,95])forthepurposesofretainingpositivityofthedensityandpressure.Wenotethatpositivity-preservinghasalsobeeninvestigatedforhydrodynamics.ThelimitersmentionedinthisparagraphcanbeviewedasgeneralizationsoflimitersthathavebeenappliedtoEuler'sequations[96,21,73].Therestofthisdissertationisasfollows.WewillbriereviewidealMHDequations7inChapter2andthenreviewtheWENOmethodologyandincludetheinterpolationandreconstructionformulasweuseforfuturereferenceinChapter3.Afterthat,inChapter4,weproposeandinvestigateourscheme,asingle-stage,single-stephighordernumericalschemeforidealMHDequations.InChapter5,weproposeandinvestigateoursecondscheme,ahighordernumericalschemeforidealMHDequations,basedonanalternativeformulationoftheWENOschemeandapplicabletocurvilinearmeshes.InChapter6,wedrawourconclusions.8Chapter2TheidealMHDequationsInthischapter,wereviewtheidealMHDequations.Inparticular,wewillreviewitshyperbolicityanddiscontinuities,bothofwhichwillbeusedlater.Intheconservationform,the(normalized)idealMHDequationsare@t26666666664ˆˆuEB37777777775+r26666666664ˆuˆuu+ptotIBBu(E+ptot)B(uB)uBBu37777777775=0;(2.1)rB=0;(2.2)whereˆisthemassdensity,ˆuisthemomentumdensity,Eisthetotalenergydensity,Bisthemagneticpisthethermalpressure,istheEuclideanvectornorm,ptot=p+12kBk2isthetotalpressure,=5=3istheidealgasconstant,andthepressuretheequationofstateE=p1+ˆkuk22+kBk22:(2.3)Herethequantitiesarenormalizedinsuchawaythatthepermeabilityandpermittivitydonotshowupexplicitlyintheequations.Inotherwords,wechooseaunitsysteminwhichtheyarebothequalto1.92.1HyperbolicityoftheidealMHDequationsWenowreviewtheofahyperbolicsystemofconservationlaws,ofwhichtheidealMHDequations(2.1)areanexample.Amorethoroughtreatmentcanbefoundin[56].WealsoreviewthewavesintheidealMHDequationsinthissection,whichcanbefoundinmostofMHDliterature.See,forexample,[48,10].2.1.Anhyperbolicsystemofconservationlawsisasystemofequationsoftheform@tq+rF=0;(2.4)wheretheconservedquantitiesq,whichisavectorfunctionoftimetandspacex,andthetensorF,whichisa(rank-2)tensorfunctionofq,aresubjecttotheconditionthatforall(unit)directionvectorn,theJacobiann@F@qisdiagonalizablewithrealeigenvalues.In2.1,ifweassumewehavencomponentsinqandmdimensionsinspaceandwriteoutthecomponents,wehaveq=(qj)j=1;:::;n;F=(Fi;j)i=1;:::;m;j=1;:::;n;n=(ni)i=1;:::;m;n@F@q=mXi=1ni@Fi;j@qkj;k=1;:::;n:(2.5)Forfuturereference,wealsowriteouttheinthendirection,whichisnF= nXi=1niFi;j!j=1;:::;m:(2.6)ThediagonalizabilityconditionoftheJacobianin2.1meansthatlocallyinanydirectionatanypointinspace,thestructureofqisapproximatelythatofplanewaves.The10eigenvaluesofn@F@qarethespeedsofthewavesindirectionnatagivenpoint.WenotehoweverthatinthecasewhenFisnotalinearfunctionofq,thesolutiontoEquation2.4canandusuallydoesdevelopverycomplexstructuresastimeevolves.IntheidealMHDsystem2.1,thewavespeedsaregivenby1;8=uncf(fastmagnetosonicwaves),(2.7)2;7=unca(Alfvenwaves),(2.8)3;6=uncs(slowmagnetosonicwaves),(2.9)4=5=un(entropyanddivergencewaves),(2.10)wherea=rpˆ(soundspeed);(2.11)ca=s(Bn)2ˆ(Alfvenspeed);(2.12)cf=0B@120B@a2+kBk2ˆ+vuut a2+kBk2ˆ!24a2(Bn)2ˆ1CA1CA12(fastmagnetosonicspeed),(2.13)cs=0B@120B@a2+kBk2ˆvuut a2+kBk2ˆ!24a2(Bn)2ˆ1CA1CA12(slowmagnetosonicspeed).(2.14)WenotethattheeigendecompositionoftheJacobianforidealMHDisverycomplicatedandisitselfasubtleissue.Theworkinthisdissertationallusesthedecompositiongiven11in[10].2.2DiscontinuitiesintheidealMHDequationsAdistinguishedfeatureofthehyperbolicsystem(2.4)isthatasolutionintheclassicalsenseusuallydonotexistforalltimet.Discontinuitiesoftendevelop,eveniftheinitialconditionsaresmooth.Therefore,oneconsidersweaksolutionstoEquation(2.4),whicharefunctionsqthatsatisfythecorrespondingintegralformofconservationlaws,andcouldpossiblycontaindiscontinuities.WenowreviewthetypesofdiscontinuitiesintheidealMHDequations(2.1).ThisknowledgewilllaterbeappliedinSection5.2.Forthehyperbolicsystemofconservationlaws(2.4),supposewestartwiththeinitialconditionsq(t=0;x)=8>>><>>>:qL;ifnx<0;qR;ifnx0;(2.15)whereqLandqRareconstantvectors,andnisanarbitrarydirectionvector.ThisiscalledaRiemannproblem.Thesolutiontosuchaproblemisafunctionqthatdependsonlyontandnx.Weareinterestedinthecasewhenthesolutionconsistsofasinglemovingdiscontinuity,thatis,whenq(t;x)=8>>><>>>:qL;ifnx<St;qR;ifnxSt;(2.16)whereSisthespeedatwhichthediscontinuitymoves.Inthiscase,theintegralformofthe12hyperbolicconservationlawmandatestheRankine-HugoniotconditionS(qRqL)=nF(qR)nF(qL):(2.17)IntheidealMHDequations,thedivergencecondition(2.2)impliesthatthemagneticintheinitialconditions(2.15)mustsatisfynBL=nBR:(2.18)Underthecondition(2.18),theRankine-Hugoniotcondition(2.17)impliesthatasinglemovingdiscontinuityinidealMHDmustbeonefromthefollowinglist.1.A(fastorslow)shock.Inthiscase,wehavenuL;nuR6=S;(2.19)nuL6=nuR;(2.20)ˆL6=ˆR:(2.21)2.Arotationaldiscontinuity.Inthiscase,wehavenuL;nuR6=S;(2.22)nuL=nuR;(2.23)ˆL=ˆR;(2.24)(uRuL)t=1pˆ(BRBL)t;(2.25)13where()trepresentsthevectorcomponenttangentialtothediscontinuityinterface(i.e.,thecomponentperpendicularton).3.Acontactdiscontinuity.Inthiscase,wehavenB6=0;(2.26)BL=BR;(2.27)uL=uR;(2.28)pL=pR:(2.29)4.Atangentialdiscontinuity.Inthiscase,wehavenB=0;(2.30)ptotL=ptotR:(2.31)HerewecanhavearbitraryjumpsinthetangentialvelocitiesandtangentialmagneticRemark2.1.Wenotethat(fastorslow)shocksaretheonlytypesofdiscontinuitiesthatcanpossiblyadmitjumpsinthenormalvelocitiesorthetotalpressures.Remark2.2.Rotational,contact,andtangentialdiscontinuitiesarelinearlydegenerate.Ro-tationaldiscontinuitiescorrespondtoAlfvenwaves.Contactandtangentialdiscontinuitiescorrespondtoentropyanddivergencewaves.Noneofthediscontinuitiesinthislistisgen-uinelynonlinear.Remark2.3.WementionedthatnBmustbethesameonthetwosidesofthediscontinu-14ity.However,inSection5.2,whenwelookatapproximateRiemannsolverstobeusedinnumericalschemesformultidimensionalproblems,wewillneedtosolveRiemannproblemswherethisconditiondoesnothold.Wewilldiscusshowwehandlethisissueinthatsection.Remark2.4.Formorediscussionsonweaksolutions,integralformsofconservationlaws,Rankine-Hugoniotconditions,lineardegeneracy,andgenuinenon-linearity,see[56].FormorediscussionsondiscontinuitiesintheidealMHDequations,see[83].15Chapter3TheWENOmethodologyInthischapter,wegiveabriefdescriptionoftheWENOinterpolationandreconstructionproceduresforfuturereferences.Thediscussionmostlyfollows[77].3.1WENOinterpolationConsiderauniformmeshwithxi=ix.Letrbeapositiveintegerandm=2r1.Supposewewanttointerpolatetheone-sidedapproximationqi+1=2fromvaluesofqonanm-point(biased)stencilS=fxir+1;:::;xi+r1g:(3.1)ApolynomialinterpolationoverSwouldgiveqi+1=2asalinearcombinationoftheval-uesfqir+1;:::;qi+r1g,whichcanthenbewrittenasaconvexlinearcombinationoftheinterpolationsobtainedfromr-pointsubstencilsofS.Tobemoreprecise,letSk=fxir+1+k;:::;xi+kg;wherek=0;:::;r1(3.2)bether-pointsubstencilsofS.ThepolynomialinterpolationoverSgivesqi+1=2=r1Xs=r+1bsqi+s=q(xi+1=2)+Ox2r1);(3.3)16wherebs'sareconstantcots.Similarly,thepolynomialinterpolationoverSkgivesq(k)i+1=2=0Xs=r+1b(k)sqi+k+s=q(xi+1=2)+Oxr);(3.4)whereb(k)s'sareconstantcots.Therethenexistconstantpositivenumbersdksuchthatqi+1=2=r1Xk=0dkq(k)i+1=2(3.5)andr1Xk=0dk=1:(3.6)IntheWENOliterature,dk'sarecalledlinearweights,becausetheapproximationqi+1=2byEquations(3.4)and(3.5)isalinearcombinationoftheqi's.Whenqisdiscon-tinuouswithinS,thepolynomialinterpolationgivesrisetoexcessiveoscillations.ThisisknownastheRunge'sphenomenonortheGibb'sphenomenon.TheWENOmethodologytacklesthisproblembyreplacingdkwithsome!kwhichr1Xk=0!k=1(3.7)and!k=dk+Oxr1);whenqissmoothonS,and!k=Oxr1);whenqisdiscontinuousonSk.(3.8)Inordertoobtain!k'ssatisfyingEquation(3.8),wecancomputesmoothnessindicators17ofthepolynomialinterpolantsoverthesubstencilsbysettingk=r1Xl=1Zxi+1=2xi1=2x2l1 @lpk(x)@lx!2dx;(3.9)wherepk(x)isthepolynomialthatinterpolatesq(k)i+1=2fromthevaluesofqonSk.Fromthesesmoothnessindicatorswecomputetheunnormalizedweightsk=dk(+k)2;(3.10)whereisasmallpositivenumber(takentobe106inallourexamples)toavoiddivisionbyzero.Normalizingthek'sgivestheWENOweights!k=kPr1s=0s:(3.11)Sincetheinterpolationgivenbyqi+1=2=r1Xk=0!kq(k)i+1=2(3.12)isnolongeralinearfunctionoftheqi's,theweights!karecallednonlinearweightsintheWENOliterature.Toillustratethealgorithmjustdescribed,let'slookattheexamplewithr=3.Inthiscase,thebigstencilisS=fxi2;xi1;xi;xi+1;xi+2g(3.13)18andthesubstencilsareS0=fxi2;xi1;xig;S1=fxi1;xi;xi+1g;S2=fxi;xi+1;xi+2g:(3.14)Equations(3.3)and(3.4)becomesqi+1=2=3128qi2532qi1+4564qi+1532qi+15128qi+2(3.15)andq(0)i+1=2=38qi254qi1+158qi;q(1)i+1=2=18qi1+34qi+38qi+1;q(2)i+1=2=38qi+34qi+118qi+2:(3.16)FromEquations(3.5),(3.15),and(3.16),wecanobtaind0=116;d1=58;d2=516:(3.17)Thesmoothnessindicators(3.9)become0=1312(qi22qi1+qi)2+14(qi24qi1+3qi)2;1=1312(qi12qi+qi+1)2+14(qi1qi+1)2;2=1312(qi2qi+1+qi+2)2+14(3qi4qi+1+qi+2)2:(3.18)Thus,theapproximationqi+1=2givenby(3.12)willsatisfyqi+1=2=q(xi+1=2)+Ox5)whenqissmoothonS,andwill,inthecasewhenScontainsadiscontinuity,beelychosentobetheinterpolationobtainedfromthesubstencilsnotcontainingthediscontinuity.19ThiscompletesthedescriptionoftheWENOinterpolationproceduretocomputeqi+1=2.Theapproximationtoq(xi+1=2)fromtheotherside,q+i+1=2,iscomputedfromthestencilS+=fxir+2;:::;xi+rg.Theinterpolationprocedureissimilarandtheformulascanbeobtainedfromtheonesforqi+1=2byexploitingthemirror-symmetrybetweenthestencilsSandS+.3.2WENOreconstructionTheWENOreconstructionsolvestheproblemofcomputingderivatives.Givenafunctionq(x),wecanaslidingfunctionhbytheimplicitformulaq(x)=1xZxx=2xx=2h(˘)d˘:(3.19)TheFundamentalTheoremofCalculusimpliesq0(x)=1xhx+x2hxx2:(3.20)Nowsupposeweknowthevaluesofqonauniformmeshxi=ix.ThisisequivalenttoknowingthevaluesoftheprimitivefunctionHofhbyH(x)=Zxx1=2h(˘)d˘;(3.21)wherex1=2isanynumber.WethuscaninterpolateHandobtainan(weighted)essentiallynon-oscillatoryinterpolantpolynomialbytheWENOinterpolationproceduredescribedinSection3.1,andtakeitsderivativetogeth(xi+1=2).Thevalueofq0(xi)can20thenbecomputedfromEquation3.20.Forstencilsofm=2r1points,theerrorinh(xi+1=2)thusobtainedisofOxm),andsoistheerrorinq0(xi),duetocancellationsintheleadingerrorterms.Inthecasewhenr=3,thebigstencilSanditssubstencilsweregivenin(3.13)and(3.14).Thevalueshreconstructedfromthesubstencilsareh(0)i+1=2=13qi276qi1+116qi;h(1)i+1=2=16qi1+56qi+13qi+1;h(2)i+1=2=13qi+56qi+116qi+2:(3.22)Thelinearweightsared0=110;d1=35;d2=310:(3.23)Thesmoothnessindicatorsk'sarethesameasin(3.18).Fromhere,wecanthencomputethenonlinearweights!k'sandtakethecombinationofthehki+1=2'susingtheseweightstoobtainhi+1=2.Forfuturereference,weshallletWENO5(qi2;qi1;qi;qi+1;qi+2)bethevalueofhi+1=2soobtained.21Chapter4Ahigh-orderpositivity-preservingsingle-stagesingle-stepmethodfortheidealmagnetohydrodynamicequationInthischapter,wedescribeourscheme,ahigh-orderpositivity-preservingsingle-stagesingle-stepschemefortheidealMHDequations.InSection4.1,wedescribetheunstaggeredconstrainedtransportframeworkweusedavectorpotentialtocontrolthedivergenceerrorofthemagneticInSection4.2,wereviewtheTaylor-discretizationPIF-WENOmethod,whichisasingle-stagemethodforsolvinghyperbolicconservationlaws.InSection4.3,wedescribehowwesolvetheevolutionequationofthevectorpotential.InSection4.4,wedescribethepositivitypreservationlimiterweused.Finally,inSection4.5,wepresentournumericalresults.Thecontentsofthischapterispublishedin[19].4.1AconstrainedtransportframeworkInthissectionwereviewtheunstaggeredconstrainedtransportframeworkproposedin[71].ThisframeworkmakesuseofthefactthatbecausethemagneticBisdivergence-free,22itadmitsavectorpotentialAthatB=rA:(4.1)Itisthereforepossibletomaintain,atthediscretelevel,thedivergence-freepropertyofBbyevolvingAalongsidetheconservedquantities,andcorrectingBtobethecurlofA.TheevolutionequationforthemagneticpotentialcanbederivedbystartingwithMaxwell'sequations.Thatis,ifwestartwithFaraday'slaw(insteadofstartingwiththeidealMHDequations),wehave@tB=E:(4.2)Next,weapproximatetheelectricbyOhm'slawforaperfectconductorE=Bu;(4.3)whichyieldstheinductionequationusedforidealMHD:@tB+r(Bu)=0:(4.4)(Afterapplyingsomeappropriatevectoridentities,Eqn.(4.4)canbewritteninaconservativeform.)Bysubstituting(4.1)into(4.4),weobtainr(@tA+(rA)u)=0:(4.5)23Thisimpliesthereisascalarfunction suchthat@tA+(rA)u= :(4.6)tchoicesof in(4.6)correspondtotgaugeconditions.Itisnotedin[43]thattheWeylgauge,whichsets 0,leadstostablenumericalsolutions.WiththeWeylgauge,theevolutionequationformagneticpotential(4.6)becomes@tA+(rA)u=0:(4.7)Intheconstrainedtransportframework,boththeconservedquantitiesq:=(ˆ;ˆu;E;B)andthemagneticpotentialAareevolved.Thetimestepfromtntotn+1proceedsasfollows:StepCT1Discretize(2.1)andupdatetheconservedquantities(ˆn;ˆun;En;Bn) (ˆn+1;ˆun+1;E;B);(4.8)whereEandBaresubjecttocorrectionsdescribedinSteps3and4below.StepCT2Discretize(4.7)andupdatethemagneticpotential:An An+1.StepCT3CorrectBtobethecurlofAthroughBn+1=rAn+1:(4.9)Notethatthisstepmodthepressure.StepCT4Modifythetotalenergysothepressureremainsunchanged(seebelow).24Inthecurrentchapter,thediscretizationweusefor(4.8)istheTaylor-discretizationofthePIF-WENOmethod[20].Thismethodisasingle-stagesingle-stepmethodforsystemsofhyperbolicconservationlaws.WegiveabriefreviewofthismethodinSection4.2.Itisnotedin[43]thatEquation(4.7)isaweaklyhyperbolicsystemthatmustbetreatedcarefullyinordertoavoidnumericalinstabilities.Alongthelinesofourpreviousworkin[23],weuseadiscretizationthatisinspiredbynumericalschemesforHamilton-Jacobiequations,inordertobeabletonon-oscillatoryderivativesofthemagneticWepresentthisdiscretizationinSection4.3.Forcertainproblems,suchasthoseinvolvinglow-plasma,negativedensityorpressureoccurinnumericalsimulations,evenwhenconstrainedtransportispresent.Thisrequirestheuseofapositivity-preservinglimiter[22,73,93].Thislimiterpreservesthepositivityofdensityandpressurebymodifyingthenumericalthatevolvestheconservedvariables.WedescribethislimiterindetailinSection4.4.Whenthislimiterisused,wealsoapplythefollowingtothepressureinStepCT4bymodifyingthetotalenergythroughEn+1=E+12(kBn+1k2kBk2):(4.10)Thisenergy\correction"keepsthepressurethesameasbeforethemagneticcorrectionin(4.9).Whilethisclearlyhasthedisadvantageofviolatingtheconservationofenergy,itisnonethelessneededinordertopreservethepositivityofpressure(otherwisethemomadetothemagnetictoobtainadivergencefreepropertycouldpotentiallycausethepressuretobecomenegative).Thistechniqueisexploredin[8,22,90].Intheresultspresentedinthecurrentwork,weusetheenergycorrection(4.10)ifandonlyifthepositivity-preservinglimiteristurnedon.254.2TheTaylor-discretizationPIF-WENOmethodWenowreviewtheTaylor-discretizationPicardintegralformulationweightedessen-tiallynon-oscillatory(PIF-WENO)method[20].Thismethodappliestogenerichyperbolicconservationlawsinarbitrarydimensions,ofwhichtheidealMHDequationisanexample.Forthepurposeofillustration,wepresentthemethodherefora2Dsystem.In2D,ahyperbolicconservationlawtakestheform@tq+@xf(q)+@yg(q)=0;(4.11)whereq(t;x;y):R+R2!Rmisthevectorofmconservedvariables,andf;g:Rm!Rmarethetwocomponentsofthefunction.Formallyintegrating(4.11)intimefromtntotn+1,onearrivesatthePicardintegralformulationof(4.11)givenbyqn+1=qnt@xFn(x;y)t@yGn(x;y);(4.12)wherethetime-averagedFandGareasFn(x;y)=1tZtn+1tnf(q(t;x;y))dt;andGn(x;y)=1tZtn+1tng(q(t;x;y))dt:(4.13)Givenadomain=[ax;bx][ay;by],thepoint-wiseapproximationsqni;jˇq(tn;xi;yj)areplacedatxi=ax+i12x;x=bxaxmx;i2f1;:::;mxg;(4.14)yi=ay+i12y;y=byaymy;y2f1;:::;myg;(4.15)26andtimet=tn,wheremxandmyarepositiveintegers.ThePIF-WENOschemesolves(4.11)bysettingqn+1i;j=qni;jtx^Fni+1=2;j^Fni1=2;jty^Gni;j+1=2^Gni;j1=2;(4.16)where^Fi1=2;jand^Gi;j1=2areobtainedbyapplyingtheclassicalWENOreconstructiontothetime-averagedFni;jandGni;jinsteadoftothe\frozen-in-time"asistradi-tionalinamethodoflines(MOL)formulation.Here,Fni;jandGni;jareapproximatedusingathird-orderTaylorexpansionintime,followedbyanapplicationofCauchy-Kovalevskayaproceduretoreplacethetemporalderivativeswithspatialderivatives.Theresultingschemeontheconservedquantitiesisthird-orderaccurateintime.Weomitthedetailshereandre-ferthereaderto[20].Extensionstothreedimensionsfollowbyincludingathirdcomponentforthefunction.4.3TheevolutionofthemagneticpotentialequationRecallthattheevolutionequationforthemagneticpotentialisin(4.7).Ifweexpandthecurlandcross-productoperator,weseethatA@tA+Nx@xA+Ny@yA+Nz@zA=0;(4.17)27whereNx=26666640uyuz0ux000ux3777775;Ny=2666664uy00ux0uz00uy3777775;Nz=2666664uz000uz0uxuy03777775(4.18)arematricesbycomponentsofthevelocityu=(ux;uy;uz).Thesystem(4.17){(4.18)hasthefollowingpropertiesthatwetakeintoaccountwhenmakingadiscretization:1.Thissystemisweaklyhyperbolic.Thismeansthatalloftheeigenvaluesof(4.17)arereal,butincontrasttobeingstrictlyhyperbolic,thesystemisnotnecessarilydiago-nalizable.Whilethisweakhyperbolicityisanartifactofhowthemagneticpotentialisevolved,itmustbetreatedcarefullytoavoidcreatingnumericalinstabilities[43];2.EachequationofthissystemcanbeviewedasaHamilton-Jacobiequationwithasourceterm[23].OneconsequenceofthissecondobservationisthatwecanuseHamilton-Jacobitechnologyinordertomagneticthathavenon-oscillatoryderivatives.Thatis,asinglederivativeofBisasecondderivativeofA,inwhichcasecareneedbetakeninordertoretainanon-oscillatoryproperty.Inthenexttwosubsections,wewilldescribehowwediscretizethesystem(4.17){(4.18),withthesetwopropertiesinmind.Beforewedothat,wemakeseveralremarks.Theisthatinthediscretization,Asitsonthesamemeshastheconservedquantities,hencethename\unstaggeredconstrainedtransport".Thesecondisthatthecurloperatorisdiscretizedinthesamewayasin[23],whichadivergencefreemagnetic(atthediscretelevel).Aremarkisthatthediscretizationwepresenthereisthird-order28accurateintime,whichisinaccordancewiththediscretizationweusefortheconservedquantities.4.3.1The2DmagneticpotentialequationWebeginwiththetwo-dimensionalcase,wheretheconservedquantitiesandthemagneticpotentialdonotdependonz.Inthiscase,thedivergence-freeconditionisrB=@xBx+@yBy=0:(4.19)ThusBzhasnoimpactonthedivergenceofB.ItthereforetocorrectonlythetwocomponentsofthemagneticAcomputationallyntwayofdoingthisistodiscardthetwocomponentsofthemagneticpotential,evolveonlythethirdcomponentaccordingtotheevolutionequation@tAz+ux@xAz+uy@yAz=0;(4.20)andforthemagneticcorrection,onlycorrectBxandBythroughBx=@yAzandBy=@xAz:(4.21)Thisisequivalenttoimposing@zAx=0and@zAy=0onthesystem(4.17){(4.18),andskippingBzinthemagneticcorrectioninEqn.(4.9).Equation(4.20)hasthefavorablepropertyofbeingstronglyhyperbolic[71],incontrasttotheweakhyperbolicityinthe3Dcase.ThisequationisalsoaHamilton-Jacobiequation,29withHamiltonprincipalfunctionAz,andHamiltonianH(t;x;y;@xAz;@yAz)=ux(t;x;y)@xAz+uy(t;x;y)@yAz:(4.22)ItisthereforepossibletodiscretizeEquation(4.20)usinganumericalschemeforHamilton-Jacobiequations[23].Inlinewiththesingle-stagesingle-stepthemeoftheproposedscheme,weuseamethodthatisinspiredbytheLax-WWENOschemesforHamilton-Jacobiequationsfrom[66].ALax-WWENOschemeforHamilton-JacobiequationsisbasedontheTaylorexpansionintimeoftheHamiltonprincipalfunction.Inordertoretainthird-orderaccuracyintime,werequireatotalofthreetimederivativesofAzi;j:Azn+1i;j=Azni;j+t@tAzni;j+t22!@2tAzni;j+t33!@3tAzni;j;(4.23)wherethetemporalderivativesontheright-hand-sidearecomputedinthefollowingmanner:1.ThederivativeisapproximatedbytheLax-FriedrichstypenumericalHamiltonian,withhigh-orderreconstructionof@xAzand@yAz[65].Namely,we@tAzi;j:=^H(@xAzi;j;@xAz+i;j;@yAzi;j;@yAz+i;j)=uxi;j @xAzi;j+@xAz+i;j2!uyi;j @yAzi;j+@yAz+i;j2!+x @xAz+i;j@xAzi;j2!+y @yAz+i;j@yAzi;j2!;(4.24)wherex=maxi;jjuxi;jjandy=maxi;jjuyi;jj(4.25)30arethemaximumvelocitiestakenovertheentiregrid,(i.e.,thenumericalistheglobalLax-Friedrichsand@xAzi;j,@yAzi;jarebyWENOreconstructionsthrough@xAzi;j=WENO5 +xAzi3;jx;+xAzi2;jx;+xAzi1;jx;+xAzi;jx;+xAzi+1;jx;!;@xAz+i;j=WENO5 +xAzi+2;jx;+xAzi+1;jx;+xAzi;jx;+xAzi1;jx;+xAzi2;jx;!;@yAzi;j=WENO5 +yAzi;j3y;+yAzi;j2y;+yAzi;j1y;+yAzi;jy;+yAzi;j+1y;!;@yAz+i;j=WENO5 +yAzi;j+2y;+yAzi;j+1y;+yAzi;jy;+yAzi;j1y;+yAzi;j2y;!;(4.26)and+xAzi;j:=Azi+1;jAzi;jand+yAzi;j:=Azi;j+1Azi;j.ThefunctionWENO5istheclassicalWENOreconstructionwhoseformulawegaveinSection3.2NotethatthisoperatorisdesignedforHamilton-Jacobiproblems,anddoesnotproducethetypicalformthatmosthyperbolicsolversproduce.2.Thehigherderivatives@2tAzand@3tAzareconvertedintospatialderivativesbywayoftheCauchy-Kovalevskayaprocedure,andtheresultingspatialderivativesareapproxi-matedusingcentralerences.Forexample,@2tAzisconvertedintospatialderivatives31by@2tAz=@t(@tAz)=@t(ux@xAzuy@yAz)=@tux@xAzux@t@xAz@tuy@yAzuy@t@yAz=@tux@xAzux@x(@tAz)@tuy@yAzuy@y(@tAz)(4.27)=@tux@xAzux@x(ux@xAzuy@yAz)@tuy@yAzuy@y(ux@xAzuy@yAz)=@tux@xAzux(@xux@xAzux@2xAz@xuy@yAzuy@x@yAz)@tuy@yAzuy(@yux@xAzux@x@yAz@yuy@yAzuy@2yAz);where@tux=@tˆuxˆ=@t(ˆux)ˆ(ˆux)@tˆˆ2;(4.28)and@tuy=@tˆuyˆ=@t(ˆuy)ˆ(ˆuy)@tˆˆ2;(4.29)with@t(ˆux),@t(ˆuy),and@tˆconvertedintospatialderivativesbywayof(2.1).NotethatwedonotitnecessarytouseaFDWENOdiscretizationforˆt.Thereareatotalof49distinctspatialderivativesthatneedtobeapproximatedin(4.27){(4.29).Similarexpressionsexistfor@3tAz.Thesearedonebyusingthecentral32formulae@xUi;j=Ui1;j+Ui+1;jx;(4.30)@2xUi;j=Ui1;j2Ui;j+Ui+1;jx2;(4.31)@x@yUi;j=Ui1;j1Ui1;j+1Ui+1;j1+Ui+1;j+1xy;(4.32)@3xUi;j=Ui2;j+2Ui1;j2Ui+1;j+Ui+2;jx3;(4.33)@2x@yUi;j=2Ui;j12Ui;j+1Ui1;j1+Ui1;j+1Ui+1;j1+Ui+1;j+1x2y;(4.34)andsimilaronesforthey-,yy-,xyy-andyyy-derivatives,whereUisanyofˆ,ˆux,ˆuy,ˆuz,E,Bx,By,Bz,orAz.Remark4.1.Inthecurrentwork,themagneticcorrectioninEqn.(4.21)isdiscretizedtofourth-orderaccuracywithcentralenceformulas.Forexample,wethey-derivativeofAztobeBxi;j=@yAzˇAzi;j28Azi;j1+8Azi;j+1Azi;j+2y;(4.35)andotherssimilarly.Thereforeinthediscretizationofequations(2.1)and(4.17),weaimforfourth-orderaccuracyinspace.Fortheequationsin(2.1),weusetheaccuratespatialdiscretizationin[20].Inthecaseof(4.17),Eqn.(4.23)isaapproximationto@tAz(tn;xi;yj)sincethereconstructions(4.26)give,forexample,@xAz+i;j=@xAz(tn;xi;yj)+Ox5);(4.36)whereasEqns.(4.30){(4.34)arechosensothatEqn.(4.23)isfourth-orderaccurateinspace33Figure4.1:Stencilneededformagneticpotential.ThevaluesofAzateach\"(savethecentermostvalue)isneededtocorrectthevalueofBatthecenter.TheshadedpointsareneededinEqns.(4.26)and(4.30){(4.34)tocomputethevalueofAzattherightmost\"point.Notethatthisstenciliscontainedinthestencilshownin[20],Fig.2.whent=Ox).Thisfourth-orderaccuracyisvinournumericalresults.Remark4.2.Ineachdimension,thestencilneededtocomputeBfromAis5pointswide,andaninspectionof(4.26)and(4.30){(4.34)indicatesthatthestencilneededtocomputeAis7pointswide.Thisresultsinastencilthatis11pointswideineachdimension.Amorecarefulanalysisof(4.30){(4.34)showsthatthisstencilisindeedcontainedwithinthestencilgivenbythePIF-WENOdiscretizationofEqn.(4.11)(whichissketchedinFig.2of[20]).WepresentasketchofthestencilrequiredtodiscretizethemagneticpotentialinFig.4.1.Wenotethatourdiscretizationsfromthatfoundin[66]inthatwedonotstoreanylowertemporalderivativesforuseincomputationofhighertemporalderivatives.Thatis,themethodproposedin[66]storesthevaluesof@tAzcomputedinEquation(4.27),andmakesuseofthecentralof@tAz.Althoughthatapproachsavesthetroubleof34expandingEqn.(4.34),ourapproachhasthefollowingadvantagesthatareinlinewithourgoals:Minimalstorage.Ourapproachavoidsthenecessityofstoringthelowertemporalderiva-tives(attheexpenseofamorecomplicatedcode).Noteherethatitwouldbenotonly@tAzand@ttAz,butalsothetemporalderivativesofˆ,ˆux,ˆuy,ˆuz,E,Bx,By,andBzthatmustbestorediftheapproachin[66]weretobetaken.Thiswouldleadtoaformidableamountoftemporarystoragewhenappliedtothemagneticpotentialevolutionequationin3D;Smallerstencils.Themethodin[66]wouldrequireastencilthatis15pointswideineachdimensioninordertoprovidethenecessaryAvaluesforcorrectingthemagneticOurmethodonlyneedsastencilthatis11pointswide.ThissomewhatthefutureworkofembeddingtheschemeintoanAMRframework,butthebiggestimprovementisthereductionoftemporarystorage.4.3.2The3DmagneticpotentialequationThe3Dmagneticpotentialequation(4.17){(4.18)isweaklyhyperbolic,andthereforespecialattentionisneeded.Thiscanbetreatedbyintroducingresistivitytermsintothesystem[23,43,44].Toillustratethetechnique,weconsidertherowofthesystemwhichis@tAxuy@xAyuz@xAz+uy@yAx+uz@zAx=0:(4.37)35Adiscretizationof@tAxinthespiritof(4.24)wouldgivethefollowingnumericallyprob-lematicformulation@tAxi;j;k=uyi;j;k @yAxi;j;k+@yAx+i;j;k2!uzi;j;k @zAxi;j;k+@zAx+i;j;k2!+y @yAx+i;j;k@yAxi;j;k2!+z @zAx+i;j;k@zAxi;j;k2!+uyi;j;k0@@xAyi;j;k+@xAy+i;j;k21A+uzi;j;k @xAzi;j;k+@xAz+i;j;k2!(4.38)wherey=maxi;j;kjuyi;j;kjandz=maxi;j;kjuzi;j;kj;(4.39)and@yAxi;j;k,@yAx+i;j;k,@zAxi;j;k,@zAx+i;j;k,@xA2i;j;k,@xA2+i;j;k,@xA3i;j;k,and@xA3+i;j;karerecon-structedinamannersimilarto(4.26).Theproblemwith(4.38)isthatthisformulationlacksnumericalresistivityinthex-direction.Asisthecasein[23],wethattheadditionofanresistivitytermto(4.38)yieldssatisfactorynumericalresults.Withtheresistivityterm,theevolutionequa-tion(4.37)becomes@tAxuy@xAyuz@xAz+uy@yAx+uz@zAx=x@2xAx;(4.40)wherexideallyx=Ox6);when@xAxissmooth;(4.41)36andx=Ox);when@xAxisnon-smooth:(4.42)Wetheresistivityxinthefollowingmanner.neasmoothnessindicatorxfor@xAxasfollows:xi;j;k=aa+a+12;(4.43)wherea=+x@xAxi;j;k)22anda+=+x@xAx+i;j;k)22:(4.44)Here,isasmallpositivenumberintroducedtoavoiddividingbyanumbercloseto0whenthepotentialissmooth(=108inallournumericalsimulations).Nowwextobex=2xx2t;(4.45)whereisapositiveconstantthatcontrolsthemagnitudeoftheresistivity.Foramoredetaileddiscussionoftheweakhyperbolicityandthereasoningleadingupto(4.45),wereferthereaderto[43,23].In[43],itisshownthathastobeintherangeof[0;0:5]tomaintainthestabilityuptoCFLonefortheirvolumemethod.Throughnumericalexperimentation,wethatinthevicinityof0:01isttocontrolpotentialoscillationsinthemagnetic37Wethususethefollowingformulationasourdiscretizationof(4.37)@tAxi;j;k=uyi;j;k @yAxi;j;k+@yAx+i;j;k2!uzi;j;k @zAxi;j;k+@zAx+i;j;k2!+y @yAx+i;j;k@yAxi;j;k2!+z @zAx+i;j;k@zAxi;j;k2!+uyi;j;k0@@xAyi;j;k+@xAy+i;j;k21A+uzi;j;k @xAzi;j;k+@xAz+i;j;k2!+2xi;j;k Axi1;j;k2Axi;j;k+Axi+1;j;kt!:(4.46)Thediscretizationof@tAyand@tAzaresimilarandweomitthemforbrevity.ForourLax-Wformulation,itremainstodiscretizethehighertemporalderivativesofthecomponentsofA.WeittoapplythesametechniquesasinSection4.3.1for@2tAzand@3tAz.Namely,weconvertallthetemporalderivativestospatialderivativesviatheCauchy-Kovalevskayaprocedure,andapproximatetheresultingspatialderivativeswithcentrals.Wenotethattheresistivitytermisonlyaddedtothetemporalderivative.4.4PositivitypreservationTheschemepresentedthusfarcanbeappliedtoalargeclassofproblems.However,forproblemswheretheplasmadensityorpressurearenearzero,Gibb'sphenomenoncancausethesevaluestobecomenegative,andhencethenumericalsimulationwillinstantlyfail.Asaingredienttothesolver,weintroduceanadditionaloptionforretainingpositivityofthesolution.Giventhelackofnumberofdegreesoffreedom,thislimitercomesattheexpenseofenergyconservation,buttheseregionsonlyoccurinsmallareaswherethedensity38orpressurebecomenegative.Apositivity-preservingschemecanbeconstructedbymodifyingthe^Fand^Gin(4.16).Let^fi+1=2;jand^gi;j+1=2bethe(global)Lax-Friedrichsby^fi+1=2;j=12(f(qni+1;j)+f(qni;j)x(qni+1;jqni;j));^gi;j+1=2=12(g(qni;j+1)+g(qni;j)y(qni;j+1qni;j));(4.47)wherexandyarethemaximalwavespeedsinthexandydirections,respectively.ThistypeofformulationiscommonlyreferredtoasLax-Friedrich'ssplittingintheENO/WENOliterature[78,79,50].TheupdateofqusingtheseLax-FriedrichsisqLFi;j=qni;jtx^fni+1=2;j^fni1=2;jty^gni;j+1=2^gni;j1=2:(4.48)Foran8-componentstatevectorq=(ˆ;ˆu;E;B),wealsointroducethenotationˆ(q)andp(q)torepresentthedensityandthethermalpressureofq.Also,letˆandpbesmallpositivenumbers.Thefollowingclaimwasconjecturedin[18].Claim4.1.Ifˆ(qni;j)>ˆandp(qni;j)>pforalli,j,andtheCFLnumberislessthanorequalto0:5,wethenhaveˆ(qLFi;j)>ˆandp(qLFi;j)>pforalli,j.Thoughnotproven,thisclaimisvin[18]usingafairlylargenumberofrandomvaluesofq.Ourpositivitylimiterassumesthisclaimistrue.However,asisnotedin[18],ifatcanbefoundsuchthatitapropertysimilartothatoftheLax-FriedrichsstatedinClaim4.1,wecanthenusethesetxesinplaceof^fi+1=2;j39and^gi;j+1=2intheconstructionofourpositivitylimiter.Themotaketheform~Fi+1=2;j=i+1=2;j(^Fi+1=2;j^fi+1=2;j)+^fi+1=2;j;(4.49)~Gi;j+1=2=i;j+1=2(^Gi;j+1=2^gi;j+1=2)+^gi;j+1=2;(4.50)wherei+1=2;jandi;j+1=2arechosensuchthat0i+1=2;j1,0i;j+1=21,theupdate(4.16)with^Fand^Greplacedbythemoes~Fand~Gin(4.49){(4.50)yieldspositivedensityandpressure,andwhilesubjecttothepositivityrequirementjuststated,i+1=2;jandi;j+1=2shouldbeascloseto1aspossible,sothatthehigh-orderareusedinregionswhereviolationofpositivityisunlikelytohappen.Following[22],wechooseeachinaseriesoftwosteps:Step-(i)Foreachi,j,\large"candidatelimitingparametersL;i;j,R;i;j,D;i;j,andU;i;j2[0;1]suchthatforall(L;R;D;U)2[0;L;i;j][0;R;i;j][0;D;i;j][0;U;i;j];(4.51)40theupdatebyq:=qni;jtx(R(^Fi+1=2;j^fi+1=2;j)+^fi+1=2;j)(L(^Fi1=2;j^fi1=2;j)+^fi1=2;j)ty(U(^Gi;j+1=2^gi;j+1=2)+^gi;j+1=2)(D(^Gi;j1=2^gi;j1=2)+^gi;j1=2)(4.52)ˆ(q)>ˆandp(q)>p;Step-(ii)Foreachi,j,seti+1=2;j:=minfR;i;j;L;i+1;jg(4.53)i;j+1=2:=minfU;i;j;D;i;j+1g:(4.54)Wenotethatwithqni;j,^F,^f,^G,and^galreadycomputed,Equation(4.52)expressestheupdateqasanfunctionofL,R,D,U,which,byabuseofnotation,isdenotedbyq(L;R;D;U).Thecotsofthe'sinq(L;R;D;U)aredenotedbyCL,CR,CD,andCU.ThuswehaveCL=tx(^Fi1=2;j^fi1=2;j);CR=tx(^Fi+1=2;j^fi+1=2;j);CD=ty(^Gi;j1=2^gi;j1=2);CU=ty(^Gi;j+1=2^gi;j+1=2);(4.55)and,with(4.48)inmind,thelimitedsolutionisq(L;R;D;U)=qLFi;j+CLL+CRR+CDD+CUU:(4.56)41Withthisnotation,Step-(i)reducestosolvingthefollowingproblemforeachi,j:Problem4.1.GivenconstantrealvectorsCL,CR,CDandCU,andconstant8-componentstatevectorqLFi;jsuchthatˆ(qLFi;j)>ˆandp(qLFi;j)>p,L;i;j;R;i;j;D;i;j;U;i;j)in[0;1][0;1][0;1][0;1]suchthatforall(L;R;D;U)in[0;L;i;j][0;R;i;j][0;D;i;j][0;U;i;j],theexpressionq(L;R;D;U)=qLFi;j+CLL+CRR+CDD+CUU(4.57)ˆ(q(L;R;D;U))>ˆ(4.58)andp(q(L;R;D;U))>p:(4.59)Theregion[0;L;i;j][0;R;i;j][0;D;i;j][0;U;i;j]shouldbe\asbigaspossible".Wenotethefollowingfact.Claim4.2.ThesetsSˆandSdbySˆ=f(L;R;D;U)2[0;1]4jˆ(q(L;R;D;U))>ˆg(4.60)42andS=f(L;R;D;U)2[0;1]4jˆ(q(L;R;D;U))>ˆandp(q(L;R;D;U))>pg(4.61)arebothconvex.Proof.Sinceˆisalinearfunctionofqandqisanefunctionof(L;R;D;U),weseethatˆ(q(L;R;D;U))isanfunctionof(L;R;D;U).ThusSˆisthepartof[0;1]4thatliesononesideofahyperplane.ThisshowsSˆisconvex.ToseetheconvexityofS,notethatbytheequationofstate2.3,thepressurepisaconcavefunctionofthecomponentsofq,wheneverˆ>0.Combinedwiththefactthatq(L;R;D;U)isanfunctionof(L;R;D;U),weseethatp(q(L;R;D;U))isaconcavefunctionof(L;R;D;U),ifˆ(q(L;R;D;U))>0.ThisshowstheconvexityofSp.NotethatProblem4.1isnotwsincethenotionof\big"isnotThealgorithmweareabouttodescribegivesasolutionthatissatisfactory,inthesensethatthisalgorithmyieldsapositivity-preservinglimiterthatbehaveswellinournumericaltests.Wenowdescribethisalgorithm.Thestepofthisalgorithmistoa\big"rectangularsubsetRˆ:=[0;ˆL][0;ˆR][0;ˆD][0;ˆU]ofSˆ.Theˆ'sarecomputed43byˆI=8>>>>>>>>>>>>><>>>>>>>>>>>>>:min8>>>>>>>>><>>>>>>>>>:1;ˆ(qLFi;j)ˆ+XJ;C(1)J<0C(1)J9>>>>>>>>>=>>>>>>>>>;ifC(1)I<0,1ifC(1)I0,(4.62)whereisasmallpositivenumber(1012inalloursimulations)andthesubscriptlettersIandJtakevaluesinL,R,D,andU.ThevalueC(1)IdenotesthecomponentoftheCIvector.ThesecondstepistoshrinkthisrectangularsubsetRˆtointoSp.TheverticesofRˆaredenotedbyAkL;kR;kD;kU,wherekI=0or1,suchthattheI-thcomponentofAkL;kR;kD;kUisAkL;kR;kD;kUI=8>>><>>>:ˆIifkI=1,0ifkI=0.(4.63)Nowforeach(kL;kR;kD;kU),shrinkAkL;kR;kD;kUtogetBkL;kR;kD;kUinthefollowingway.Ifp(AkL;kR;kD;kU)p,setBkL;kR;kD;kU=AkL;kR;kD;kU.Otherwise,wesolveforthesmallestpositiversuchthatp(rAkL;kR;kD;kU)p,andsetBkL;kR;kD;kU=rAkL;kR;kD;kU.Inordertosolveforthisvalue,weapplyatotalof10iterationsofthebisectionmethod.Amoret(approximate)solvercouldeasilyreplacethisstep(e.g.,asinglestepofthemethodoffalseposition).NotethatthetverticesAkL;kR;kD;kUareingeneralshrunkbytfactorsr.NowwesetI;i;j=min(kL;kR;kU;kD);kI=1BkL;kR;kU;kDI;(4.64)44wherethesubscriptIindicatestheI-thcomponent.NotethatthisisequivalenttoarectangularsubsetinsidetheconvexpolygonwithverticesBkL;kR;kD;kU,kI=0;1.Thiscompletesthedescriptionofourpositivity-preservinglimiterin2D.The3Dcaseissimilar.4.5NumericalresultsInthissection,wepresenttheresultsofnumericalsimulationsusingtheproposedmethod.Unlessotherwisespconstrainedtransportisturnedon,thegasconstant=5=3,theCFLnumberis0:5,andfor3Dsimulationstheviscositycotis=0:01.4.5.1SmoothAlfvenwaveThesmoothAlfvenwaveproblemisoftenusedforconvergencestudiesofnumericalschemesforidealMHDequations[43,71,90].Thisproblemisaone-dimensionalproblem(computedinmultipledimensions)thathasaknownsmoothsolution.In1D,theinitialconditionsforthisproblemare(ˆ;ux;uy;uz;uz;p;Bx;By;Bz)(0;x)=(1;0;0:1sin(2ˇx);0:1cos(2ˇx);0:1;1;0:1sin(2ˇx);0:1cos(2ˇx)):(4.65)Theexactsolutionto(4.65)propagateswiththeAlfvenspeedthatisunity(i.e.,q(t;x)=q(0;x+t)).The2Dand3DsmoothAlfvenwaveproblemsareobtainedfromthe1Dproblembyrotatingthedirectionofwavepropagation.454.5.1.1SmoothAlfvenwave:The2DproblemThe2DversionofthesmoothAlfvenwaveproblemisobtainedbyrotatingthedirec-tionofpropagationbyanangleof˚,sothatthewavenowpropagatesindirectionn=cos˚;sin˚;0i.Identicalto[43,23],thecomputationaldomainweuseis[0;1=cos˚][0;1=sin˚],where˚=tan1(0:5).Periodicboundaryconditionsareappliedonallfoursides.Forthisproblem,wepresentnumericalresultswithandwithouttheenergycorrection.InTable4.1,weobservetheoverallthird-orderaccuracyofthemethod,andinTable4.2,wetfasterthanthemeshspacingaswellasrunthesolutiontoashortertimeinordertoextractthespatialorderofaccuracy.Forthesetestcases,weobservethepredictedfourth-orderaccuracyinspace.ofconvergenceinspace,third-orderintime,andlittlebetweentheresultsobtainedwithandwithouttheenergycorrectionturnedon.Whentheforthepositivity-preservinglimiteristurnedoninthecode,weseeidenticalresultsaswithoutit,becausethisproblemdoesnothavedensityorpressurethatisnearzero.Thechoiceofmy=2mxallowsforx=y.4.5.1.2SmoothAlfvenwave:The3DproblemThesetupweusehereisthesameasthatusedin[43],Section6.2.1.Thedirectionofpropagationisn=cos˚cos;sin˚cos;sini;(4.66)andthecomputationaldomainis0;1cos˚cos0;1sin˚cos0;1sin;(4.67)46MeshCFLErrorinBOrderErrorinAzOrder32640.53:842105|5:356106|641280.54:9401062:967:5301072:831282560.56:3241072:979:6971082:962565120.58:0201082:981:2181082:99(a)MeshCFLErrorinBOrderErrorinAzOrder32640.53:848105|5:320106|641280.54:9381062:967:4691072:831282560.56:3181072:979:6281082:962565120.58:0091082:981:2101082:99(b)Table4.1:2DsmoothAlfvenwave.Here,weshowL1-errorsatatimeoft=1:0forthesolutionwithandwithouttheenergy\correction"step.Table(a)hastheenergycorrectionturnedandTable(b)hastheenergycorrectionturnedon.Becausetimeisonlydiscretizedtothird-orderaccuracy,weobservethepredictedthird-orderaccuracyofthesolverhere.MeshCFLErrorinBOrderErrorinAzOrder32640.53:852106|1:078107|641280.252:3561074:038:1211093:731282560.1251:4661084:015:19010103:972565120.06259:11710104:013:29110113:98(a)MeshCFLErrorinBOrderErrorinAzOrder32640.53:852106|1:078107|641280.252:3561074:038:1211093:731282560.1251:4661084:015:19010103:972565120.06259:11710104:013:29110113:98(b)Table4.2:2DsmoothAlfvenwave.Here,weshowL1-errorsatashorttimeoft=0:01forthesolutionwithandwithouttheenergy\correction"step.Table(a)hastheenergycorrectionturnedandTable(b)hastheenergycorrectionturnedon.Here,wetfasterthanxinordertoexposethespatialorderofaccuracyofthesolver.Becauseweonlyuseafourth-orderaccuratespatialdiscretizationforrA,weonlyobservefourth-orderaccuracydespitethefactthatthevariablesarediscretizedtoaccuracy.47MeshCFLErrorinBOrderErrorinAOrder1632320.54:784104|5:116105|3264640.52:4521054:293:1811064:01641281280.53:0931062:994:6121072:791282562560.53:9691072:966:1331082:91(a)MeshCFLErrorinBOrderErrorinAOrder1632320.54:882104|5:176105|3264640.52:4851054:303:1911064:02641281280.53:1051063:004:6211072:791282562560.53:9771072:966:1471082:91(b)Table4.3:3DsmoothAlfvenwave.InthistableweshowtheL1-errorsatamoderatetimeoft=1:0.InTable(a)thepositivity-preservinglimiterisandthethepositivity-preservinglimiter(andtheenergycorrectionstep)isturnedonfortheresultsinTable(b).Becausetimeisdiscretizedtothird-orderaccuracy,themethodisformallyonlythird-orderaccurateintime.InTable4.4werunthesolvertoashorttimeinordertoexposethespatialorderofaccuracy.where˚==tan1(0:5).Periodicboundaryconditionsareimposedonalldirections.Weagainseektonumericallyinvestigatethespatialandtemporalordersofaccuracy,withandwithouttheenergycorrectionstep.TheerrorsinBandAarepresentedinTables4.3{4.4.Herewechoosemy=mz=2mxsothatx=y=z=cos.Similartothe2Dcase,weobservefourth-orderofconvergenceinspace,third-orderintime,andlittlebetweentheresultsobtainedwithandwithouttheenergycorrectionstepturnedon.4.5.22DrotatedshocktubeproblemSimilartothesmoothAlfvenproblems,therotatedshocktubeproblemisa1Dproblem,withdirectionofwavepropagationrotatedacertainangle.Thesetupweuseinthecurrent48MeshCFLErrorinBOrderErrorinAOrder1632320.56:752105|5:715107|3264640.254:2801063:983:8561083:89641281280.1252:6661074:002:6131093:881282562560.06251:6521084:011:71110103:93(a)MeshCFLErrorinBOrderErrorinAOrder1632320.56:752105|5:715107|3264640.254:2801063:983:8561083:89641281280.1252:6661074:002:6131093:881282562560.06251:6521084:011:71210103:93(b)Table4.4:3DsmoothAlfvenwave.Here,weshowtheL1-errorsatashorttimeoft=0:01.Inaddition,wetfasterthanthemeshspacinginordertoextractthespatialorderofaccuracy.Table(a)hasthepositivity-preservinglimiterturnedandTable(b)hasthepositivity-preservinglimiter(aswellasthecorrectionstep)turnedon.Theresultsarealmostidentical.workisthesameasthatin[15],whichwerepeathereforcompleteness.Theinitialconditionsconsistofashock(ˆ;u?;uk;uz;p;B?;Bk;Bz)=8>>><>>>:(1;0;0;0;1;0:75;1;0)if˘<0,(0:125;0;0;0;0:1;0:75;1;0)if˘0,(4.68)where˘=xcos˚+ysin˚,andu?andB?arevectorcomponentsperpendiculartotheshockinterface,andukandBkarethevectorcomponentsparalleltotheshockinterface.Namelyux=u?cos˚uksin˚;uy=u?sin˚+ukcos˚;(4.69)Bx=B?cos˚Bksin˚;By=B?sin˚+Bkcos˚:(4.70)49TheinitialconditionformagneticpotentialisAz(0;x;y)=8>>><>>>:0:75˘if˘0,0:75+˘if˘>0,(4.71)where=xsin˚+ycos˚.Thecomputationaldomainis[1;1][0:5;0:5]witha200100mesh.Theboundaryconditionsusedarezerothorderextrapolationontheconservedquantitiesandextrapolationonthemagneticpotential.Thatis,wesettheconservedquantitiesattheghostpointstobeidenticaltothelastvalueontheinteriorofthedomain,andwevaluesforthemagneticpotentialatghostpointsthroughrepeatedextrapolationoftwopointstencilsstartingwithtwointeriorpoints.Onthetopandbottomboundaries,thedirectionofextrapolationisparalleltotheshockinterface.InFigure4.2wepresentresultsforthedensityofsolutionscomputedusingPIF-WENOwithandwithoutconstrained-transport.WenotethatthecontourplotofthesolutionobtainedwithoutconstrainedtransportdoesnotexhibittheunphysicalwigglesasisthecaseinFig2(b)of[23].However,ascanbeseenfromtheplotsofthesliceaty=0,unphysicaloscillationsappearinthePIF-WENOschemethathasconstrainedtransportturnedItisalsoclearfromtheseplotsthattheconstrainedtransportmethodweproposeinthecurrentworkisabletosuppresstheunphysicaloscillationssatisfactorily.Asasidenote,weithelpstouseaglobal,asopposedtoalocalvalueforintheLax-Friedrichsxsplittingforthehigh-orderWENOreconstructioninordertofurtherreduceundesirablespuriousoscillations.50Figure4.2:Therotatedshocktubeproblem.Here,wecomparethesolverwith(leftpanels)andwithout(rightpanels)constrainedtransportturnedon.Auniformmeshofsize200100isusedforbothsimulations.Theangleofrotationfortheinitialconditionsis˚=tan1(0:5),and30equallyspacedcontoursrangingfromtheminumumtothemaximumofeachfunctionareusedforthetoptwopanels.ThecontourplotforthesolutionwithoutCTcontainssmallwigglesthataremuchmorepronouncedwhenslicesofthesolutionaresampled.Tothisend,asliceofthesolutionalongy=0ispresentedinthebottomtwopanels.FurtherevolutionproducesasolutionthatcausesthecodetofailinthecasewhereCTisturned.Thepositivity-preservinglimiteristurnedinordertoexercisethecode.Thesolidlinesinthebottomimagesarereferencesolutionsthatarecomputedbysolvingtheequivalent1Dshockproblemonauniformmeshwith50,000pointswiththeWENOmethod.51WithConstrainedTransportNoConstrainedTransportFigure4.3:2Drotatedshocktubeproblem.Theleftpanelshaveconstrainedtransportturnedon,andtherightpanelshaveconstrainedtransportturnedComponentsofthemagneticatt=0:2alongtheslicey=0.Themeshsizeis200100.Thepositivity-preservinglimiteristurnedinordertoexercisethecode.EachsolidlineisthereferencesolutiondescribedinFig.4.2.Weobservethatconstrainedtransport(withtheHamilton-Jacobisolver)allowsustonumericallycomputemagneticwithfarfeweroscillationsthanwouldotherwisebeobtainable.524.5.32DOrszag-TangvortexInthissection,weinvestigatetheOrszag-Tangvortexproblem.Anotablefeatureofthisproblemisthatshocksandcorticesemergefromsmoothinitialconditionsastimeevolves.ThisisastandardtestproblemfornumericalschemesforMHDequations[23,26,71,90,94].Theinitialconditionsareˆ=2;u=(sin(y);sin(x);0);B=(sin(y);sin(2x);0);p=;(4.72)withaninitialmagneticpotentialofAz(0;x;y)=0:5cos(2x)+cos(y):(4.73)Thecomputationaldomainis[0;2ˇ][0;2ˇ],withdouble-periodicboundaryconditions.WepresentinFigure4.4thedensitycontourplotsatt=0:5,t=2,t=3,andt=4,computedbyPIF-WENOwithconstrainedtransportonandpositivity-preservinglimiteron.Controlofdivergenceerrorofthemagneticeldiscriticalforthistestproblem.Ifweturntheconstrainedtransport,thesimulationcrashesatt=1:67.InFigure4.5wepresentplotsofthedensityattimet=1:5obtainedwithtofthenumericalschemes.Notethedevelopmentofthenonphysicalfeaturesinthesolutionobtainedwithoutconstrainedtransport.Wealsonotethatinthisproblem,wherethepositivity-preservinglimiterisnotneededforthesimulation,thelimiterleadstolittleinthesolution.WepresentinFigures4.5and4.6plotsthatdemonstratethis.TheplotsofthetypeinFigures4.4and4.6are53(a)(b)(c)(d)Figure4.4:TheOrszag-Tangproblem.Weshowdensitycontourplotsat(a)t=0:5,(b)t=2,(c)t=3,(d)t=4.ThesolutioniscomputedusingPIF-WENOwithconstrainedtransportonandpositivity-preservinglimiteron,ona192192mesh.Atotalof15equallyspacedcontoursareusedforeachgraph.54presentedinmanysources.Ourplotsagreewithresultsfromtheliterature[23,26,71,90,94].4.5.42DrotorproblemThisisatwo-dimensionaltestproblemthatinvolveslowpressurevalues[8].Thesetupweusehereisthesameastheverylowversionfoundin[91].Theseinitialconditionsareˆ=8>>>>>>>><>>>>>>>>:10ifr0:1,1+9~f(r)ifr2(0:1;0:115),1ifr0:115,(4.74)ux=8>>>>>>>><>>>>>>>>:10y+5ifr0:1,(10y+5)~f(r)ifr2(0:1;0:115),0ifr0:115,(4.75)uy=8>>>>>>>><>>>>>>>>:10x5ifr0:1,(10x5)~f(r)ifr2(0:1;0:115),0ifr0:115,(4.76)uz=0;Bx=2:5p4ˇ;By=0;Bz=0;p=108;Az=2:5p4ˇy;(4.77)wherer=q(x0:5)2+(y0:5)2;~f(r)=13(23200r):(4.78)Thecomputationdomainweuseis[0;1][0;1],withzerothorderextrapolationontheconservedquantitiesandorderextrapolationonthemagneticpotentialastheboundary55(a)(b)(c)(d)Figure4.5:TheOrszag-Tangproblem.Densityplotsatt=1:5forsolutionscomputedusingtinthenumericalscheme.(a)PIF-WENO,withconstrainedtransportturned(b)PIF-WENOwithconstrainedtransportandpositivitypreservinglimiterturnedon;(c)PIF-WENOwithconstrainedtransportturnedonandpositivitypreservinglimiterturned(d)Theslicealongy=0:5072.Thesolutionsarecomputedwitha192192mesh.Atotalof15equallyspacedcontoursareusedforeachof(a),(b),and(c).56(a)(b)(c)Figure4.6:TheOrszag-Tangproblem.Pressureplotsatt=3:0forsolutionscomputedusingtinthenumericalscheme.(a)PIF-WENOwithconstrainedtransportturnedonandpositivitypreservinglimiterturned(b)PIF-WENOwithconstrainedtransportandpositivitypreservinglimiterturnedon;(c)Theslicealongy=1:9799.Thesolutionsarecomputedwitha192192mesh.Atotalof15equallyspacedcontoursareusedforeachof(a)and(b).57(a)(b)Figure4.7:The2Drotorproblem.Shownhereare(a)thepseudocolorplotofthemagneticpressureand(b)themagneticline.Thissolutioniscomputedonameshofsize400400.Constrainedtransportandpositivity-preservinglimitersareturnedon.conditionsonallfoursides(i.e.,conservedquantitiesattheghostpointsaresetequaltothelastinteriorpoint,andvaluesforthemagneticpotentialareedthroughrepeatedextrapolationoftwopointstencils).Wecomputethesolutiontoatimeoft=0:27usinga400400meshandpresenttheplotsofthemagneticpressureandthemagneticeldlineinFigure4.7.ThemagneticlineplotpresentedhereisthecontourplotofAz.Atotalof50levelsareusedinthiscontourplot.Ourresultisconsistentwiththeonepresentedin[91].Wenotethatthepositivity-preservinglimiterisnecessarytocompletethistestproblem,becauseotherwisethecodefails.584.5.5Cloud-shockinteractionThecloud-shockinteractionproblemisastandardtestproblemforMHD[23,26,43,44,71].Theinitialconditionsare(ˆ;ux;uy;uz;p;Bx;By;Bz)=8>>>>>>>><>>>>>>>>:(3:86859;11:2536;0;0;167:345;0;2:1826182;2:1826182)ifx<0:05;(10;0;0;0;1;0;0:56418958;0:56418958)ifx>0:05andr<0:15,(1;0;0;0;1;0;0:56418958;0:56418958)otherwise,wherer=q(x0:25)2+(y0:5)2in2D,andr=q(x0:25)2+(y0:5)2+(z0:5)2in3Ddenotesthedistancetothecenterofthestationarycloud.Forbothproblems,weusetheinitialmagneticpotentialAx=0;Ay=0;Az=8>>><>>>:2:1826182x+0:080921431ifx0:05;0:56418958xifx0:05:(4.79)Inthe2Dcase,weonlykeeptrackofAz.4.5.5.1Cloud-shockinteraction:The2DproblemThecomputationaldomainweuseis[0;1][0;1],withzerothorderextrapolationontheconservedquantitiesandorderextrapolationonthemagneticpotentialastheboundaryconditionsonallfoursides(i.e.,conservedquantitiesattheghostpointsaresetequaltothelastinteriorpoint,andvaluesforthemagneticpotentialareedthroughrepeatedextrapolationoftwopointstencils).59(a)(b)(c)(d)(e)(f)Figure4.8:The2Dcloud-shockinteractionproblem.Here,werunthesolvertoatimeoft=0:06.Inthethreepanels,weshowSchlierenplotsfor(a)thenaturallogofthedensity,(b)thenormofthemagneticand(c)thepressureforameshofsize256256.Thesameresultsforameshofsize512512arepresentedinpanels(d){(f),whereweobservemuchhigherresolutionfortheproblem.Constrainedtransportandpositivity-preservinglimiterareturnedonforbothsimulations.60Wecomputethesolutionatt=0:06usinga256256mesh.TheSchlierenplotsoflnˆandofjBjarepresentedinFigure4.8.Wenoteherethatthecurrentschemeisabletocapturetheshock-wave-likestructurenearx=0:75.Thisisconsistentwithourpreviousresultin[23],andanimprovementoverearlierresultsin[26,71].Wealsonotethatthepositivity-preservinglimiterisnotrequiredforthissimulation.Nonetheless,wepresenttheresultheretodemonstratethehighresolutionofourmethod,evenwhenthelimiteristurnedon.4.5.5.2Cloud-shockinteraction:The3DproblemThecomputationaldomainforthisproblemis[0;1][0;1][0;1],withzerothorderextrap-olationontheconservedquantitiesandorderextrapolationonthemagneticpotentialastheboundaryconditionsonallsixfaces(i.e.,conservedquantitiesattheghostpointsaresetequaltothelastinteriorpoint,andvaluesforthemagneticpotentialarethroughrepeatedextrapolationoftwopointstencils).Wecomputethesolutiontoatimeoftimet=0:06ona256256256mesh.InFigure4.9,weshowtheevolutionofthedensityofthesolution.Wehavetworemarksonthisresult.Theisthattheshock-wave-likestructurenearx=0:75atthetimeisalsovisiblewhenweusea128128128mesh.Thesecondisthatthepositivity-preservinglimiterisrequiredtorunthissimulationwith256256256mesh.Thereasonisthatextrastructurethatcontainsverylowpressureshowsupinthismeshattimet=0:0378.Thisextrastructurecannotbeobservedonthecoarsermeshwiththemethodproposedinthecurrentwork,nordoweobserveitwithourpreviousSSP-RKsolver[23].Therefore,itdoesnotcausetroubleforsimulationsusinga128128128mesh.61Figure4.9:The3Dcloud-shockinteractionproblem.Schlierenplotsofln(ˆ).Thesolutionhereiscomputedusinga256256256mesh.Cross-sectionsaty=0:5andz=0:5fortheregiony0:5andz0:5areshown.Constrainedtransportandpositivity-preservinglimiterareturnedon.624.5.6BlastwaveexampleIntheblastwaveproblems,strongshocksinteractwithalow-background,whichcancausenegativepressureifnothandledproperly.Theseproblemsareoftenusedtotestthepositivity-preservingcapabilitiesofnumericalmethodsforMHD[5,8,22,34,57,63,97].Theinitialconditionscontainapiecewisepressure:p=8>>><>>>:0:1r<0:1;1000otherwise;(4.80)whereristhedistancetotheorigin,andaconstantdensity,velocityandmagnetic(ˆ;ux;uy;uz;Bx;By;Bz)=(1;0;0;0;100=p4ˇ=p2;100=p4ˇ=p2;0):(4.81)TheinitialmagneticpotentialissimplyA=(0;0;100y=p4ˇ=p2100x=p4ˇ=p2):(4.82)In2DweonlykeeptrackofAz,aswedowithallthe2Dexamples.4.5.6.1Blastwaveexample:The2DproblemInthissection,wepresentourresultonthe2Dversionoftheblastwaveproblem.Thecomputationaldomainis[0:5;0:5][0:5;0:5],withzerothorderextrapolationontheconservedquantitiesandorderextrapolationonthemagneticpotentialastheboundaryconditionsonallfoursides(i.e.,conservedquantitiesattheghostpointsaresetequaltothelastinteriorpoint,andvaluesforthemagneticpotentialareedthroughrepeated63extrapolationoftwopointstencils).Resultsforthesolutioncomputedtoatimeoft=0:01ona256256mesharepresentedinFigure4.10.There,wedisplaycontourplotsofˆ,p,juj,andjBj.Theseplotsarecomparabletothepreviousresultsin[22].Wenotethatnegativepressureoccursrightinthestepifpositivity-preservinglimiteristurned4.5.6.2Blastwaveexample:The3DproblemForthe3Dversionoftheblastwaveproblem,wechoosethecomputationaldomaintobe[0:5;0:5][0:5;0:5][0:5;0:5]withzerothorderextrapolationontheconservedquantitiesandorderextrapolationonthemagneticpotentialastheboundaryconditionsonallsixfaces(i.e.,conservedquantitiesattheghostpointsaresetequaltothelastinteriorpoint,andvaluesforthemagneticpotentialarethroughrepeatedextrapolationoftwopointstencils).Thesolutionatt=0:01iscomputedusinga150150150mesh.WepresentinFigure4.11theplotsofthedensityandpressure,andalsoinFigure4.12thecontourplotsofthesliceatz=0ofthedensity,pressure,velocity,andmagneticpressure.Theseresultsarecomparabletothosefoundin[22,34,63,97].Wenoteherethatnegativepressureoccursinthesecondtimestepifthepositivity-preservinglimiteristurned4.5.7ErrorsinenergyconservationWhenthepositivity-preservinglimiteristurnedon,wemakeuseofanenergycorrectionstepinEqn.(4.10)inordertokeepthepressurethesameasbeforethemagneticcorrection.Thisbreakstheconservationoftheenergy,andthereforeweinvestigatetheectofthisstepforseveraltestproblems.Because(global)energyconservationonlyholdsforproblemsthat64(a)(b)(c)(d)Figure4.10:2Dblastproblem.Shownherearethecontourplotsatt=0:01of(a)density,(b)thermalpressure,(c)magnitudeofvelocity,and(d)magneticpressure.Atotalof40equallyspacedcontoursrangingfromthemintothemaxofthefunctionareusedforeachplot.Themeshsizeis256256.65(a)(b)Figure4.11:3Dblastproblem.Shownherearethepseudocolorplotsatt=0:01of(a)densityand(b)pressure.Themeshsizeis150150150.InFigure4.12weplotacutofthesolutionalongz=0.Thepositivity-preservinglimiterisrequiredtosimulatethisproblem.haveeitherperiodicboundaryconditionsorconstantvaluesneartheboundarythroughouttheentiresimulation,weonlychooseproblemswiththispropertyforourtestcases.Forthe2Dproblems,the(relative)energyconservationerrorattimet=tnisasEnergyconservationerror:=Pi;jEni;jE0i;jPi;jE0i;j;(4.83)andtheenergyconservationerrorsforthe3Dproblemsarenedsimilarly.Resultsforthe2Dand3DsmoothAlfventestcasearepresentedinFigure4.13,whereweobservenegligibleerrorsproducedbytheenergycorrectionstep.Weattributethistothefactthatthisproblemretainsasmoothsolutionfortheentiretyofthesimulation.ResultsforproblemswithshocksandvorticesarepresentedinFigure4.14,wherewenon-zeroerrors.Foreachofthesetestproblems,wepresenttheresultsfromseveraltsizesofmeshes.AlltheproblemsareruntothetimefoundinSections4.5.1{4.5.6saveone.66(a)(b)(c)(d)Figure4.12:3Dblastproblem.Shownarethecontourplotsattimet=0:01cutatz=0of(a)density,(b)thermalpressure,(c)normofvelocity,and(d)magneticpressure.Thesolutionisobtainedusinga150150150mesh.Atotalof40equallyspacedcontoursareusedforeachplot.67(a)(b)Figure4.13:EnergyconservationerrorsforthesmoothAlfventestcases.Shownhereareresultsforthe(a)2DsmoothAlfven,and(b)3DsmoothAlfventestcases.Inordertoextracttheerrors,weplottheresultsonasemi-logscalebecauseotherwisetheresultsareindiscerniblefromthet-axis.Forthissmoothtestcase,theofthetheenergycorrectionstep(andhencethepositivity-preservinglimiter)isnegligiblebecausethesolutionremainssmoothfortheentiresimulation.Forthe2DOrszag-Tangproblem,werunthesimulationstoamuchlatertimeoft=30inordertoquantifytheenergyconservationerrorsforalongtimesimulationonanon-trivialproblem.Finally,inFigure4.15wealsoincludetheconservationerrorswhenthepositivity-preservinglimiteristurnedWeobservethatthesolverretainstotalenergyuptomachineerrors,asshouldbethecase.Wenotethefollowingpatternsintheenergyconservationerrors:Theerrorsarebelow1%forallthetestproblemsinthedurationofthesimulationspresented;Whenthepositivity-preservinglimiteristurnedon,theerrorsgrowlinearlyintimeanddecreaseasthemeshisWethereforeconcludethattheviolationinenergyconservationintroducedbythepositivity-preservinglimiteristfortheproblemstestedinthiswork.68(a)(b)(c)(d)Figure4.14:Energyconservationerrors.Shownhereareconservationerrorswhenthepositivity-preservinglimiter(andhencetheenergycorrectionstep)isturnedon.(a)2DOrszag-Tangproblem,(b)2Drotorproblem,(c)2Dblastproblem,and(d)3Dblastproblem.Notethattherotorandblastproblemsrequiretheapplicationofapositivity-preservinglimiterinordertorun,butthiscomesattheexpenseoflosingenergyconservation.FortheOrszag-Tangtestproblem,weruntoalatetime.Again,weobservethattheerrorsinenergyconservationdecreaseasthemeshis(a)(b)Figure4.15:Energyconservationerrors.Here,weshowresultsforconservationerrorswhenthepositivity-preservinglimiter(andhencetheenergycorrectionstep)isturnedThesolveranalyticallyconservesthediscretetotalenergyuptomachineprecision.Shownhereareresultsforthe(a)2DsmoothAlfven,and(b)3DsmoothAlfventestcases.Notethelogarithmicscalefortheaxes,andthattheseerrorsarenumericallyconserveduptomachineprecision.69Chapter5Ahigh-orderschemeforidealmagnetohydrodynamicsoncurvilinearmeshesbasedonanalternativeformulationoftheWENOschemeandconstrainedtransportInthischapter,wedescribeoursecondscheme,ahigh-ordererenceschemefortheidealMHDequations,whichisbasedonanalternativeformulationoftheWENOschemeandusesconstrainedtransporttocontrolthedivergenceerrorofthemagneticInSection5.1,wedescribeanalternativeformulationoftheWENOscheme,whichisbasedontheschemeproposedin[51,52],withtheadditionofanlimiteronthehigherorderterms.ThisschemewillbethebaseschemeweusetosolvetheidealMHDequations.InSection5.2,wedescribeseveralHLLtypeapproximateRiemannsolversfortheidealMHDequations.Inparticular,wedescribetheHLLDsolverin5.2.3,whichisoneofthetwoRiemannsolversweexperimentinacomponentofourbasescheme(seeEquation(5.22)),theotherbeingLax-Friedrichs.InSections5.3and5.4,wedescribeconstrainedtransportand70positivity-preservinginthecurrentcontext.FinallyinSection5.5,wepresentthenumericalresults.5.1AnalternativeformulationoftheWENOschemeInthissection,wedescribeaWENOschemebasedonanalternativeformulationof5.1.1IntroductionThesimplesthyperbolicconservationlawtakestheform@q@t+@f(q)@x=0;(5.1)whereqisascalarfunctionoftandx,andfisafunctionofq.Hereqistheconservedquantityandfisthe.NotethisisEquation(2.4)withm=n=1.InordertosolveEquation(5.1),weuseasemi-discreteconservativeerencescheme,whereweuseauniformmeshwithxi=ixanddiscretizeEquation(5.1)inspacetoget@qi@t+1x(^fi+1=2^fi1=2)=0:(5.2)Hereqisitsonxiand^fissomenumericalthatsitsonhalfgridpointsanddependsontheqi's.Thisnumerical1x^fi+1=2^fi1=2=@xf(q(x))jxi+Oxm);(5.3)71wheremisthespatialorderofaccuracyofthescheme.Equation(5.1)isthenintegratedintimeusingsomeODEsolver,e.g.Runge-Kuttamethod.Inthecurrentwork,wehaveusedthethird-orderSSP-RKmethodproposedin[39].Thenumerical^fisahighorderapproximationtotheso-calledslidingfunctionhoff(q).Ifwelethbethefunctionsuchthatf(q(x))=1xZxx=2xx=2h(˘)d˘;(5.4)wewouldhave^fi+1=2=h(xi+1=2)+Oxm):(5.5)Equations(5.4)and(5.5)actuallysaysthatwearegiventheaverageofhovercellscenteredatthegridpointsxi'sandweneedtoapproximatethevaluesofthefunctionhathalfgridpointsxi+1=2'sfromthosecellaverages.Thisisthereconstructionproblemattheheartofstandardvolumeschemes[41,60]andhasbeenusedtoconstructvariousnumericalschemesforhyperbolicconservationlaws[79,49,6,77].Itispossibletocompute^fi+1=2directlyfromthevaluesoffongridpointsusingtheproceduredescribedinSection3.2.ThisissometimecalledtheWENOscheme.Wereferthereaderto[77]fordetails.Itwasnotedin[76]thatwhileEquation(5.3)requirestheerrorinEquation(5.5)beofOxm+1),theOxm)erroristbecausewhenfissmooth,cancellationsintheerrortermswillgiveriseto(5.3).Itisnotedin[78]thataTaylorexpansionoftheslidingfunctionhatxi+1=2givesthe72followingexpressionofhi+1=2intermsofthecellaveragefunctionf,hi+1=2=fi+1=2+[(m1)=2]Xk=1a2kx2k @2k@x2kf!i+1=2+Oxm);(5.6)wherethea2k'saresomeconstants.Intheexamplesweshowlater,weuseatruncationatm=5,andthereforeapproximate^fby^fi+1=2=fi+1=2124x2@2xfji+1=2+75760x4@4xfji+1=2:(5.7)TheterminEquation(5.7)isapproximatedbyfi+1=2=F(qi+1=2;q+i+1=2);(5.8)whereFisamonotone[86]andqi+1=2aretlyhigh-orderone-sidedapproxima-tionstoqatxi+1=2.TheWENOinterpolationdescribedinSection3.1isusedtoobtainqi+1=2.ThehigherordertermsinEquation(5.7)areapproximatedusingcentraltogetherwithalimiterconstructedfromtheWENOsmoothnessindicatorsin(3.18).Inwhatfollows,weshalldescribetheaforementionedlimiterinSection5.1.2.Afterthat,wewilldescribehowwehandleone-dimensionalsystems,multidimensionalcases,andcurvilinearmeshesinSections5.1.3,5.1.4,and5.1.5,respectively.5.1.2AlimiteronthehigherordertermsinthenumericalInthissection,welookatapproximationstothetermsinthesummationinEquation(5.6).InorderthatweretainanOxm)error,theapproximationto@2k@x2kfi+1=2mustbe73oferrorOxm2k).Inthespecialcasewhenm=5,thisrequirementcanbebyusingthecentralx2@2xfji+1=2ˇ148(5fi2+39fi134fi34fi+1+39fi+25fi+3);x4@4xfji+1=2ˇ12(fi23fi1+2fi+2fi+13fi+2+fi+3):(5.9)NotethatthestencilneededinEquation(5.9)iscontainedintheunionofSusedintheWENOinterpolationproceduredescribedinSection3.1.Equation(5.9)istheapproachusedin[52]toapproximatethesecondandthirdtermsin(5.7).Intheoriginalarticle[51]onthealternativeformulationWENOscheme,thehigherorderderivatives@2xfji+1=2and@4xfji+1=2arefurtherexpandedintermsofthederivativesoffwithrespecttoqandthederivativesofqwithrespecttox,beforecentralareapplied.Inbothofthesearticles,whichdealwithEulerequationsofhy-drodynamics,satisfactoryresultscanbeobtainedwithoutusingspecialtreatmentsontheseterms.However,inexperimentingwiththeMHDequations,weitnecessarytoapplyanadditionallimiterontheseterms.Inthecaseofm=5,wewouldliketomultiplythesecondandthirdtermsontherighthandsideof(5.9)byanumber˙that˙=1+Ox3);whenqissmoothonthestencilS=fxi2;:::;xi+3g;˙=Ox2);whenqcontainsastrongdiscontinuityonS.(5.10)Such˙canbeconstructedfromthesmoothnessindicatorsusedinWENOinterpolation74inthefollowingway.Westartfromthek'sinEquation(3.18)andthenset˙max=1+j02j+minf0;1;2g;˙min=1+j02j+maxf0;1;2g;(5.11)whereisasmallpositivenumber(takentobe106inallourexamples)toavoiddivisionbyzero.Wecanthusobtainacandidateforthecot˙by˙=˙min˙max:(5.12)Notethat˙dependsonlyonthestencilS=fxi2;:::;xi+2g.AsimilarformulaforS+=fxi1;:::;xi+3ggivesrisetoanothercandidate˙+.Wethenset˙=minf˙;˙+g:(5.13)5.1.3One-dimensionalsystemsWeconsiderone-dimensionalhyperbolicsystemsofconservationlaws,whicharethespecialcaseof(2.4)withm=1.Aone-dimensionalsystemofconservationlawtakestheform@q@t+@f(q)@x=0;(5.14)whereq=(q1(t;x);:::;qn(t;x))(5.15)75isavectorfunctionoftandxandf(q)=(f1(q);:::;fn(q))(5.16)isthefunction.Heresuperscriptsdenotethecomponentsofthevectors.Supposethesystem(5.14)ishyperbolic,thatis,theJacobian@f=@qhasnrealeigenvalues1(q)n(q)(5.17)andasetofnindependent(right)eigenvectorsr1(q);:::;rn(q):(5.18)Thus,ifweletR(q)=(r1(q);:::;rn(q))(5.19)bethematrixwhosecolumnsaretheeigenvectorsof@f=@q,wewillhaveR1(q)@f@qR(q)=diag(1(q);:::;n(q)):(5.20)SimilartoEquation(5.7),wehave,inthecaseofspatialorderofaccuracym=5,thenumerical^fi+1=2=fi+1=2124x2@2xfji+1=2+75760x4@4xfji+1=2;(5.21)76wherefi+1=2=F(qi+1=2;q+i+1=2)(5.22)isobtainedfromsomemonotoneuxFresultedfromanyofavarietyofRiemannsolvers,withqi+1=2obtainedfromWENOinterpolation,andthehigherordertermsisobtainedusingacombinationofcentralandanadditionallimiter,sameasdescribedinSection5.1.2.Inthecaseofahyperbolicsystem,theWENOinterpolationweusetoobtainqi+1=2isperformedonthelocalcharacteristicvariablesinsteadofonthecomponentsofq.Tobemoreprecise,wehavethefollowingalgorithm.1.Ateachhalfgridpointxi+1=2,computeanaveragestateqi+1=2.InalltheexamplesweshowlaterfortheMHDequations,wehaveusedthearithmeticmeanofprimitivevariables,thatis,set i+1=2=12( i+ i+1)(5.23)forthevariables rangingoverfˆ;u;p;Bg,andrecovertheconservedvariablesq=fˆ;ˆu;E;Bgfromthe 's.2.Ateachhalfgridpointxi+1=2,computetherightandlefteigenvectorsoftheJacobian@f=@q.3.Ateachhalfgridpointxi+1=2,projecttheconservedquantitiesinthestencilneededforcomputingthenumerical^fi+1=2ontothelocalcharacteristicvariables.Thatis,letvi;k=R1i+1=2qi;fork=r+1;:::;r;(5.24)whereweusevi;ktodenotethelocalcharacteristicvariables,andR1istheleft77eigenvectorscomputedinlaststep.4.Foreachi,performtheWENOinterpolationonthevi;k'stogetvi+1=2.5.Foreachi,projectvi+1=2backontotheconservedquantities,qi+1=2=Ri+1=2vi+1=2:(5.25)5.1.4MultidimensionalcasesHyperbolicconservationlawsinmorethanonespatialdimensionscanbetreatedinadimension-by-dimensionfashion,asisusuallydoneinschemes.Forexample,asystemofhyperbolicconservationlawintwodimensionstakestheform@q@t+@f(q)@x+@g(q)@y=0;(5.26)whereqisavectorfunctionoft,x,andy,andfandgaretheinthexandydirections,respectively.Onauniformmeshwithxi=ixandyj=jy,Equation(5.26)canbesolvedbyusingthesemi-discretescheme@qi;j@t+1x(^fi+1=2;j^fi1=2;j)+1y(^gi;j+1=2^gi;j1=2)=0;(5.27)where1x^fi+1=2;j^fi1=2;j=@xf(q(x;y))j(xi;yj)+Oxm1)(5.28)and1y^gi;j+1=2^gi;j1=2=@yg(q(x;y))j(xi;yj)+Oym1)(5.29)78Figure5.1:Aschematicdiagramofthephysicalandcomputationaldomainforcurvilinearmeshes.foraschemewithspatialorderofaccuracyequaltom.Now^fand^gcanbeapproximatedbyusingthemethoddescribedintheprevioussections.5.1.5CurvilinearmeshesInthissection,weconsiderthesituationwherethemeshonthephysicaldomainisobtainedfromauniformmeshonthecomputationaldomainviaacontinuousmap.Supposewehaveatwo-dimensionalscalarequationofhyperbolicconservationlaw(5.26).Also,supposethecoordinates(x;y)isrelatedtothecurvilinearcoordinates(˘;)viasomecontinuouscoordinatetransformation.WealsoletJ=˘x˘yxy=x˘xy˘y1(5.30)betheJacobianofthetransformation.SeeFigure5.1foraschematicdiagram.79TransformingEquation(5.26)intoanequationofconservationlawinthecurvilinearcoordinatesgives@eq@t+@ef@˘+@eg@=0;(5.31)whereeq=q=J;ef=e˘xf+e˘yg;eg=fxf+eyg;(5.32)forthemetrictermsbye˘x=˘x=J=y;e˘y=˘y=J=x;fx=x=J=y˘;ey=y=J=x˘:(5.33)WenotethatEquation(5.31)ishyperbolicaslongasEquation(5.26)ishyperbolic.Ifweequipthedomaininthe(˘;)-planewithauniformmesh,wecansolveEqua-tion(5.31)usingthemethodwedescribedintheprevioussections.Inthecourseofdoingthat,wewillneedtocomputethemetrictermsinEquation(5.33),whichcanbeapproxi-matedusingntlyhighordercentral5.2HLLtypeRiemannsolversforMHDequationsInthissection,wedescribeseveralHLLtypeRiemannsolversforMHDequations,whichweshalluseinthealternativeformulationoftheWENOscheme.SeeEquation(5.22).TheHLLapproximateRiemannsolverwasproposedbyHarten,Lax,andvanLeer[42].ThissolveremploysanapproximatesolutiontotheRiemannprobleminwhichthereisoneintermediatestateconnectedtotheleftandrightstatesbydiscontinuities.Thisinterme-diatestateisobtainedbyexploitingtheconservationoftheequations,whichiscalledthe80consistencyconditionintheliterature.TheHLLsolverappliestoanyhyperbolicconserva-tionlawanddoesnotmakeuseofknowledgeontheparticularsystem.WhenappliedtotheEulerequationsforhydrodynamics,theHLLsolverexhibitsexcessivedissipationinthepresenceofcontactdiscontinuities.Toremedythis,Toro,Spruce,andSpeares[88]proposedwhatisnowknownastheHLLC(Cforcontact)solverfortheEuler'sequations.Thissolverassumestwointermediatestatesintheapproximatesolution,connectedtoeachotherbyacontactdiscontinuityandtotheleftandrightstatesbyshocks.TheRankine-Hugoniotcondition,inadditiontotheconsistencycondition,isusedtodeterminetheintermediatestates.SimilarideaswerealsoappliedtotheidealMHDequationsin[40,58,64].Gurski[40]andLi[58]namedtheirsolvers\HLLC"solversforMHD,becausetheyassumetwointerme-diatestatesconnectedtoeachotherbyacontactdiscontinuity,aswasthecaseintheHLLCsolversfortheEuler'sequations.MiyoshiandKusano[64]namedtheirsolver\HLLD"(Dfordiscontinuities)solver,becausetheirsolverassumesfourintermediatestatesandcanexactlyresolvemosttypesofdiscontinuitiesintheidealMHDequations,theonlyexceptionbeingtheslowshocks.TheHLLDsolverisoneofthetwoRiemannsolversweexperimentwhenweexploretheofthechoiceofRiemannsolverinthealternativeformulationoftheWENOscheme,theotherbeingLax-Friedrichs.In[64],itisassumedthatthenormalvectortothediscontinuityinterfaceisparalleltooneofthecoordinateaxesandthatthemagneticcomponentnormaltothediscontinuityinterfaceisthesameatallpointsinspace.NeithernecessarilyholdstrueintheproblemsweneedtheRiemannsolvertohandle.CurvilinearmeshesdemandthatweconsiderRiemannproblemswithinitialconditions(2.15)wherethenormalvectorntodiscontinuityinterfacecouldbeanydirection.Multidimen-sionalMHDproblemsallowsforjumpsinthenormalmagneticldbecausedivergence-free81conditioncanbesavedbysuchjumpsinotherdirections.WepresentourversionofanHLLDsolverthatcanhandletheseissuesbydroppingtheassumptionsonthedirectionnandthenormalmagneticSincetheHLLDsolverusescertainresultsobtainedfromtheHLLandHLLCsolvers,wepresentthesetwosimplersolversinSections5.2.1and5.2.2.Afterthat,wepresenttheHLLDsolverinSection5.2.3.5.2.1TheHLLapproximateRiemannsolverConsidertheRiemannproblemgivenbytheinitialconditions(2.15).WeuseFtodenotetheinthendirectionwhichwewroteasnFin(2.6).Wealsointroducethenotation[a]=a+a(5.34)todenotethejumpofaquantityaacrossadiscontinuityinterfacenormalton,wherea+isthevalueonthesidethatnpointsto,andaisonthesidethatnpointsfrom.Forexample,forthediscontinuityintheinitialconditions(2.15),[ˆ]=ˆRˆL.Theapproximatesolution~qemployedintheHLLsolverconsistsofthreestates,qR,qHLL,andqR,separatedbytwodiscontinuitiespropagatingatwavespeedsSLandSR,respectively,whereSLSR.Inotherwords,~q(t;xn)=8>>>>>>>><>>>>>>>>:qL;if(xn)=tSL;qHLL;ifSL(xn)=tSR;qR;ifSR(xn)=t:(5.35)HereSLandSRarerespectivelythesmallestandthelargestofallthesignalspeeds,ap-82proximatedbythesmallestandlargestamongallwavespeedsinqLandqR.Thein-termediatestateqHLLissothatforallT,˘L,and˘RsuchthatT>0and˘Lmin(0;TSL)max(0;TSR)˘R,theconsistencyconditionZ˘R˘L~q(t;˘)d˘=˘RqR˘LqL+TFLTFR(5.36)is(see[86]).Inthesubsonicsituation,whereSL<0<SR,theconsistencyconditionwith˘L=TSLand˘R=TSRimpliesqHLL=SRqRSLqL+FLFRSRSL:(5.37)TheGodunovtypenumericalcorrespondingtothisapproximatesolutionisF=8>>>>>>>><>>>>>>>>:FL;if0SL;SRSL+SLSR(qRqL)SRSL;ifSL0SR;FR;ifSR0:(5.38)5.2.2TheHLLCapproximateRiemannsolverTheHLLCapproximateRiemannsolverwepresentinthissectionwasproposedbyLiin[58].Theapproximatesolutionemployedinthissolverhastwointermediatestates83connectedbyacontactdiscontinuity.Thatis,theapproximatesolutionisgivenby~q(t;xn)=8>>>>>>>>>>>><>>>>>>>>>>>>:qL;if(xn)=tSL;qL;ifSL(xn)=tSM;qR;ifSM(xn)=tSR;qR;ifSR(xn)=t;(5.39)whereSLandSRaretheestimatedsmallestandlargestsignalspeeds,SMistheestimatedspeedoftheentropywave(aswellasthatofthedivergencewave),andqLandqRaresomeintermediatestatestobedetermined.Oncethisapproximatesolutionisdetermined,thecorrespondingGodunovtypenumericalisF=8>>>>>>>>>>>><>>>>>>>>>>>>:FL;if0SL;FL=FL+SL(qLqL);ifSL0SM;FR=FR+SR(qRqR);ifSM0SR;FR;ifSR0:(5.40)WenowproceedtodescribetheapproximatesolutioninthecasewhereSL0SR.WementionedbeforethatthenormalvelocitiesandthetotalpressuresdonotchangeacrossacontactdiscontinuityinMHD.ItisthereforereasonabletoassumeuLn=uRn=SM=uHLLn(5.41)84andptotL=ptotR:(5.42)WeshallfurthermaketheassumptionthatBLn=BRn=BHLLn:(5.43)TheRankine-Hugoniotcondition(2.17)acrossS,whereiseitherLorR,impliesSqF=SqF:(5.44)Equation(5.44)onˆgivesSˆˆun=Sˆˆun:(5.45)Thereforeˆ=ˆSunSSM;(5.46)inviewof(5.41).Equation(5.44)onˆugivesSˆu(ˆ(un)u+ptotn(Bn)B)=Sˆu(ˆ(un)u+ptotn(Bn)B):(5.47)85Takingdotproductsofbothsidesof(5.47)withn,wegetSˆSMˆS2Mptot+(Bn)2=SˆSMˆ(un)2ptot+(Bn)2:(5.48)Fromthisand(5.46),wegetptot=ptot+ˆ(Sun)(SMun)+(Bn)2(Bn)2:(5.49)Thusweareabletocomputefrom(5.48)themomentaintheintermediatestates,ˆu=ˆ(Sun)+(ptotptot)n+(Bn)B(Bn)BSSM(5.50)ToobtainthetotalenergyEintheintermediatestates,weagainusetheRankine-Hugoniotcondition(5.44)acrossS.ThistimewegetSE((E+ptot)(un)(uB)(Bn))=SE((E+ptot)(un)(uB)(Bn)):(5.51)SolvingforE,wegetE(Sun)+ptotSMptot(un)+(Bn)(Bu)(Bn)(Bu)SSM:(5.52)WestillneedtodetermineBin(5.50)and(5.52).Weshallusetheconsistencycondition,whichinthecurrentcontexttakestheformSMSLSRSLqL+SRSMSRSLqR=qHLL:(5.53)86Apply(5.53)and(5.37)toˆuandsimplifywith(5.42)takenintoaccount,andwegetBRnBR+BLnBL=0:(5.54)Inviewof(5.43),itisreasonabletosetBL=BR=BHLL:(5.55)Nowtheintermediatestatesintheapproximatesolutionarecompletelydeterminedin(5.41),(5.43),(5.46),(5.49),(5.55),(5.50),and(5.52).Li[58]notedthatiftheRankine-Hugoniotcondition(5.44)onBwereusedinplaceoftheconsistencycondition(5.53)onˆutogetB,theconsistencyconditionwillbeviolatedandinstabilitieswillhappen.5.2.3TheHLLDapproximateRiemannsolverTheapproximatesolutionemployedintheHLLDsolverhasfourintermediatestates,con-nectedbytworotationaldiscontinuitiesandonecontactortangentialdiscontinuity.Inotherwords,~q(t;xn)=8>>>>>>>>>>>>>>>>>>>>><>>>>>>>>>>>>>>>>>>>>>:qL;if(xn)=tSL;qL;ifSL(xn)=tSL;qL;ifSL(xn)=tSM;qR;ifSM(xn)=tSR;qR;ifSR(xn)=tSR;qR;ifSR(xn)=t;(5.56)87whereSLandSRareestimatedfastshockspeeds,SLandSRareestimatedspeedsoftherotationaldiscontinuities,andSMistheestimatedspeedofthecontactortangentialdiscontinuity.ThecorrespondingGodunovtypenumericalisF=8>>>>>>>>>>>>>>>>>>>>><>>>>>>>>>>>>>>>>>>>>>:FL;if0SL;FL=FL+SLqLSLqL;ifSL0SL;FL=FL+SLqLSLSLqLSLqL;ifSL0SM;FR=FR+SRqRSRSRqRSRqR;ifSM0SR;FR=FR+SRqRSRqR;ifSR0SR;FR;ifSR0:(5.57)WenowproceedtodescribetheapproximatesolutioninthecasewhereSL0SR.Similartotheassumptions(5.41),(5.42),and(5.43)madeinHLLCsolver,weassumeuLn=uLn=uRn=uRn=SM=uHLLn;(5.58)ptotL=ptotL=ptotR=ptotR;(5.59)andBLn=BLn=BRn=BRn=BHLLn:(5.60)RecallthatrotationaldiscontinuitiesarelinearlydegenerateandcorrespondtotheAlfvenwaves,whichpropagateatspeedsunq(Bn)2=ˆ.Since[ˆ]=0acrossrotational88discontinuities,wecansetSL=SMs(BHLLn)2ˆL(5.61)andSR=SM+s(BHLLn)2ˆR;(5.62)whereˆcanbeobtainedbyusingtheRankine-Hugoniotcondition(5.44)onˆacrossS,whichgivesˆ=ˆSunSSM;(5.63)thesameas(5.46).Thefour-intermediate-stateapproximatesolution(5.56)degeneratestoatwo-intermediate-statesolutionwhenSLandSRareeithertooclosetoSLandSR,respectively,orwhentheyaretooclosetoSM.Inbothcases,wefallbacktotheHLLCsolverandforegowhatwedescribebelow.Wenowconsiderthenon-degeneratecasewhereSL,SL,SL,SM,SR,SR,andSRaretlyspreadapart.ThestateqisobtainedfromtheRankine-HugoniotconditionacrossS.Thevalueofˆwasgivenin(5.63).Equations(5.49)and(5.50)inHLLCremainvalidinthecurrentcontext.Werepeatthemhere:ptot=ptot+ˆ(Sun)(SMun)+(Bn)2(Bn)2;(5.64)ˆu=ˆ(Sun)+(ptotptot)n+(Bn)B(Bn)BSSM:(5.65)89NowapplytheRankine-HugoniotconditionacrossSonB,whichgivesSB((un)B(Bn)u)=SB((un)B(Bn)u):(5.66)Withungivenby(5.58),Bngivenby(5.60),ˆgivenby(5.63),andptotgivenby(5.64),weseethat(5.65)and(5.66)formalinearsystemof(algebraic)equationsinuandB.Thesolutiontothissystemisu=ˆ(Sun)(SSM)(BHLLn)(Bn)ˆ(Sun)(SSM)(BHLLn)2u+(ptotptot)(SSM)ˆ(Sun)(SSM)(BHLLn)2n+(Bn)(SSM)(BHLLn)(Sun)ˆ(Sun)(SSM)(BHLLn)2B(5.67)andB=ˆ(Sun)2(BHLLn)(Bn)ˆ(Sun)(SSM)(BHLLn)2B+ˆ(Sun)(BnBHLLn)ˆ(Sun)(SSM)(BHLLn)2u(ptotptot)(BHLLn)ˆ(Sun)(SSM)(BHLLn)2n:(5.68)UsingtheRankine-HugoniotconditionacrossSonE,wecanobtainE=E(Sun)+ptotSMptot(un)+(Bn)(Bu)(Bn)(Bu)SSM;(5.69)whichhasthesameformas(5.52).WenotethatwhereasintheHLLCsolverweusedconsistencycondition(5.53)onˆutoobtainaequationforB(5.55),intheHLLDsolver90weusedtheRankine-HugoniotconditionacrossSonˆuandBtogetˆuandB,andsaveaconsistencyconditionsimilarto(5.53)fordetermining(ˆu)andBinthetwoinnerintermediatestates.Wenowproceedtothevaluesoftheinnerintermediatestatesq.TheRankine-HugoniotconditionacrossSonˆgivesSˆ(ˆu)n=Sˆ(ˆu)n:(5.70)Solve(5.70)with(5.58)takenintoaccount,wegetˆ=ˆ:(5.71)Notethatthisisconsistentwiththefactthat[ˆ]=0acrossarotationaldiscontinuity,whichweassumedinderivingtheestimatesofSLandSRin(5.61)and(5.62).SincewehavebeenconsideringthecasewhereSL,SL,SL,SM,SR,SR,andSRaretlyspreadapart,wecanassumethatBHLLn6=0;(5.72)inviewof(5.61)and(5.62).Undertheassumption(5.72),theRankine-HugoniotconditionacrossSMonˆuimpliesBL=BR(5.73)91andthesameconditiononBimpliesuL=uR:(5.74)Theconsistencycondition(consistencywiththeintegralconservationlaw)inthecurrentcontexttakestheformSRSRqR+SRSMqR+SMSLqL+SLSLqLSRqR+SLqL+FRFL=0:(5.75)Weregroupthetermsin(5.75)andgetSRSMqR+SMSLqL+SRqRqR+FRSLqLqL+FLSRqR+SLqL=0:(5.76)TheRankine-Hugoniotcondition(2.17)acrossSgivesS(qq)=FF(5.77)Substitute(5.61),(5.62),(5.77)into(5.76),andwegetjBHLLnj0@qRqˆR+qLqˆL1A+FRFLSRqR+SLqL=0:(5.78)Applying(5.78)toˆuwith(5.71)and(5.74)takenintoaccount,wegetjBHLLnjqˆR+qˆLu+((un)ˆu+ptotn(Bn)B)R((un)ˆu+ptotn(Bn)B)LSR(ˆu)R+SL(ˆu)L=0:(5.79)92Solvingforuwith(5.58),(5.59),(5.60),(5.61),and(5.62)takenintoaccountgivesu=qˆLuL+qˆRuR+sign(BHLLn)BRBLqˆL+qˆR:(5.80)Applying(5.78)toBwith(5.73)takenintoaccountgivesjBHLLnj0@1qˆR+1qˆL1AB+((un)B(Bn)u)R((un)B(Bn)u)LSRBR+SLBL=0:(5.81)SolvingforBwith(5.58),(5.60),(5.61),and(5.62)takenintoaccountgivesB=qˆRBL+qˆLBR+sign(BHLLn)qˆLˆRuRuLqˆL+qˆR:(5.82)TheRankine-HugoniotconditionacrossSonEgivesSE((E+ptot)(un)(uB)(Bn))=SE((E+ptot)(un)(uB)(Bn)):(5.83)SolvingforEwith(5.58),(5.59),(5.60),(5.61),and(5.62)takenintoaccountgivesE=Epˆ(uBuB)sign(BHLLn);(5.84)whereand+correspondsto=Land=R,respectively.Nowtheintermediatestatesintheapproximatesolutionarecompletelydeterminedin(5.58),(5.59),(5.60),(5.61),(5.62),(5.63),(5.64),(5.67),(5.68),(5.69),(5.71),(5.73),93(5.74),(5.80),(5.82),and(5.84).5.3ConstrainedtransportInourcurrentscheme,weusethesameconstrainedtransportframeworkasinSection4.1tocontrolthedivergenceerrorofthemagneticNamely,weevolveavectorpotentialA,which(4.1)B=rA:Theevolutionequationisgivenby(4.7)@tA+(rA)u=0:Inthecurrentmultistagesetting,wesetB=rAaftereachstageinthetimeintegrator.AsexplainedinSection4.3.1,inthecaseoftwospatialdimensions,thedivergence-freeconditionbecomes(4.19)rB=@xBx+@yBy=0;whereBzdoesnotplayanyrole.ItthereforecestocorrectonlyBxandByin2D,whichcanbedonebyevolvingAzaccordingto(4.20)@tAz+ux@xAz+uy@yAz=094andcorrectingBxandBythrough(4.21)Bx=@yAz;By=@xAz:Inthecurrentmultistagescheme,Equation(4.20)issolvedbyaWENOmethodforHamilton-Jacobiequationsinthesamewayasin[23].Themagneticcorrection(4.21)isdiscretizedbyusingfourth-orderaccuracycentralnces,alsointhesamewayasin[23].Constrainedtransportoncurvilinearmeshesisdonesimilarly.Inthecurvilinearcoordi-nates,Equation(4.20)becomes@tAz+J(uxe˘x+uye˘y)@˘Az+J(uxfx+uyey)@Az=0(5.85)andEquation~(4.21)becomesBx=J(e˘y@˘Az+ey@Az);By=J(e˘x@˘Az+fx@Az):(5.86)Equation(5.85)issolvedbyaWENOmethodforHamilton-JacobiequationsandEqua-tion(5.86)isdiscretizedbyusingfourth-orderaccuracycentralWenotethatwhilesuchdiscretizationonlyguaranteesthedivergence-freeconditionofmagnetictotruncationerrors,inpracticewethisisttosuppresstheunphysicaloscillationsassociatedwiththedivergenceerrorofB.955.4Apositivity-preservinglimiterWhenweapplythecurrentschemewiththeHLLDRiemannsolver,weitnecessarytoapplyapositivity-preservinglimitersimilartotheonedescribedinSection4.4.ThisispossiblyduetotheenhancedresolutionprovidedbytheHLLDsolver.WeencounteredasimilarsituationinSection4.5.5.2,whereameshrevealsaregionoflowpressurenotresolvedincoarsermeshes.Thepositivity-preservinglimitergoesmostlythesameasthelimiterinSection4.4.Theareasfollows.1.Onacurvilinearmesh,theequationactuallybeingsolvedisgivenby(5.31).There-foretheexpressionsfortheupdatesandinSection4.4needtobemoaccordingly.2.Inthecurrentmultistagescheme,theproceduredescribedinSection4.4istobeappliedateachstage.Wenotethatwhileitissuggestedin[22]thatitisttoapplythelimiteronlyatthestageofeachtimestep,weitnecessarytoapplyitateachstageinthecurrentscheme.This,again,ispossiblyduetotheenhancedresolutionprovidedbytheHLLDsolver.3.Wehaveusedfandgtodenotethehighorderinthecurrentchapter.InChapter4,however,weusedFandGtodenotethehigh-ordertime-averagedandfandgtodenotethelow-orderItshouldbeclearhowtotranslatethedescriptioninSection4.4intoadescriptionwiththenotationsinthecurrentchapter.965.5Numericalresults5.5.12DsmoothAlfvenwaveproblemWeusethesmoothAlfvenwaveproblemtotesttheorderofconvergenceofournumericalscheme.Theinitialconditionsweregivenin(4.65)as(ˆ;ux;uy;uz;uz;p;Bx;By;Bz)(0;x)=(1;0;0:1sin(2ˇx);0:1cos(2ˇx);0:1;1;0:1sin(2ˇx);0:1cos(2ˇx)):TheexactsolutionpropagateswiththeAlfvenspeed,whichequals1(i.e.,q(t;x)=q(0;x+t)).InSection4.5.1werotatedthedirectionofAlfvenwavepropagationsothatitisnotparalleltoanymeshlines.Inthecurrentsection,wekeepthedirectionofwavepropagationparalleltothex-axisandtweakthemeshtoachievethesameThecomputationaldomainissettobe(˘;)2[0;1]2,withthemeshlinesperturbedaccordingtox=˘+xsin(2ˇax);y=+ysin(2ˇ˘ay);(5.87)wherexandyarethemagnitudeofperturbationandaxandayarethewavenumbersoftheperturbation.Intheresultsweshowbelow,weusex=0:01,y=0:02,ax=2,anday=4.Asanillustration,weplota3232meshwiththeseparametersinFigure5.2.InTables5.1and5.2,wepresenttheL1-errorsinBandAz,obtainedusingthealter-nativeWENOmethodwiththeLax-FriedrichsandwiththeHLLDrespectively.Theseresultsthatourmethodhasfourthorderofaccuracy.Wenotethatforthis97Figure5.2:Perturbedmeshusedfor2DAlfvenproblem.3232mesh.(˘;)2[0;1]2.x=˘+0:01sin(2ˇ2),y=+0:02sin(2ˇ˘4).MeshErrorinBOrderErrorinAzOrder32325:131103|1:560104|64644:0901043.658:3941064.221281282:6581053.945:2661073.992562561:6771063.993:3221083.99Table5.1:L1-errorsin2DsmoothAlfvenproblem.AlternativeWENOwithLax-FriedrichsConstrainedtransportandpositivity-preservinglimiteron.smoothproblem,thebetweenusingtheLax-FriedrichsandtheHLLDisverysmall.5.5.2Brio-WushocktubeThisisafrequentlytestedRiemannproblem.Theinitialconditionsareasgivenin(4.68),(ˆ;ux;uy;uz;p;Bx;By;Bz)=8>>><>>>:(1;0;0;0;1;0:75;1;0)ifx<0,(0:125;0;0;0;0:1;0:75;1;0)ifx0.(5.88)98MeshErrorinBOrderErrorinAzOrder32325:257103|1:604104|64644:0901043.688:3881064.261281282:6581053.945:2641073.992562561:6771063.993:3221083.99Table5.2:L1-errorsin2DsmoothAlfvenproblem.AlternativeWENOwithHLLDux.Constrainedtransportandpositivity-preservinglimiteron.InSection5.5.2.1,werunourschemeonbothuniformandnon-uniform1Dmeshes,soastotestitshandlingofshocksacrossabruptchangesinmeshspacing.InSection5.5.2.2,werunourschemeona2Dmesh,withtheinitialconditionsrotated,soastotestthehandlingofshocksandcontrolofdivergenceerrorofthemagneticld.5.5.2.11DshocktubeInFigures5.3{5.6,wepresentplotsofthedensityandthetransversemagneticusingvarioussetupsof1Dsimulations.Inallresults,wehaveusedmeshesof200points.ThemeshesinFigures5.3and5.4areuniform.ThemeshesinFigures5.5and5.6aredeterminedbyx=8>>><>>>:109˘ifj˘j0:2,sign(˘)29+209(j˘j0:2)otherwise,(5.89)with1˘1.WenotethatboththealternativeWENOmethodswithLax-FriedrichsandwiththeHLLDcanhandletheabruptchangeinthespacinginthenon-uniformcaseswithoutintroducingtoomuchoscillations.99Density.TransversemagneticFigure5.3:Brio-Wushocktube.t=0:2.AlternativeWENOwithLax-FriedrichsAuniformmeshof200pointswereused.Positivity-preservinglimiteristurnedon.Density.TransversemagneticFigure5.4:Brio-Wushocktube.t=0:2.AlternativeWENOwithHLLDx.Auniformmeshof200pointswereused.Positivity-preservinglimiteristurnedon.100Density.TransversemagneticFigure5.5:Brio-Wushocktube.t=0:2.AlternativeWENOwithLax-FriedrichsAnon-uniformmeshof200pointswereused.Positivity-preservinglimiteristurnedon.Density.TransversemagneticFigure5.6:Brio-Wushocktube.t=0:2.AlternativeWENOwithHLLDAnon-uniformmeshof200pointswereused.Positivity-preservinglimiteristurnedon.1015.5.2.22DrotatedshocktubeWealsoruntheBrio-Wushocktubeteston2Dmeshes.Inalltheresultspresentedinthissection,weusea200100uniformmeshonthedomain[1;1][0:5;0:5].Wehaverotatedtheinitialconditionbyanangleoftan1(0:5).Inthissetup,theenforcingofthedivergence-freeconditionbecomescrucial.Withouttheconstrainedtransport,thesimulationusingalternativeWENOwithHLLDcrashesbeforetreaches0:2,whereasthesimulationusingalternativeWENOwithLax-Friedrichsdemonstratesexcessiveunphysicaloscillations,whichcanbeseenfromtheplotsin5.7.Ontheotherhand,theplotsinFigures5.8and5.9showthatturningonconstrainedtransportcansatisfactorilycontroltheunphysicaloscillations.5.5.32DOrszag-TangvortexInthissection,wetestourschemeonthe2DOrszag-Tangvortexproblem,whichweinves-tigatedinSection4.5.3withourscheme.Theinitialconditionsweregivenin(4.72),ˆ=2;u=(sin(y);sin(x);0);B=(sin(y);sin(2x);0);p=;(5.90)withaninitialmagneticpotentialgivenin(4.73),Az(0;x;y)=0:5cos(2x)+cos(y):(5.91)102Contourplotsofdensity.Densityalongy=0.Transversemagneticalongy=0.Figure5.7:RotatedBrio-Wushocktube.t=0:2.AlternativeWENOwithLax-Friedrichs200100uniformmesh.Constrainedtransportandpositivity-preservinglimiterturned103Contourplotsofdensity.Densityalongy=0.Transversemagneticalongy=0.Figure5.8:RotatedBrio-Wushocktube.t=0:2.AlternativeWENOwithLax-Friedrichs200100uniformmesh.Constrainedtransportandpositivity-preservinglimiterturnedon.104Contourplotsofdensity.Densityalongy=0.Transversemagneticalongy=0.Figure5.9:RotatedBrio-Wushocktube.t=0:2.AlternativeWENOwithHLLD200100uniformmesh.Constrainedtransportandpositivity-preservinglimiterturnedon.105Figure5.10:Orszag-Tangtestproblem.Contourplotofdensityatt=3.Aperturbed192192mesh.15equallyspacedcontourlines.AlternativeWENOwithLax-FriedrichsConstrainedtransportandpositivity-preservinglimiterturnedon.Thecomputationaldomainis(˘;)2[0;2ˇ][0;2ˇ],andthemeshisgivenbyx=˘+xsin(ax);y=+ysin(˘ay);(5.92)wherex=0:03,y=0:05,ax=2,anday=4.Wepresentthecontourplotsofthedensityattimet=3inFigures5.10and5.11.Theresultsagreewiththosefoundintheliterature[23,26,71,90,94].Wenotethatwehavealsosuccessfullyrunthesimulationtothemuchlatertimeoft=10,whichshowstherobustnessoftheourscheme.106Figure5.11:Orszag-Tangtestproblem.Plotofdensityatt=3.Aperturbed192192mesh.15equallyspacedcontourlines.AlternativeWENOwithHLLDConstrainedtransportandpositivity-preservinglimiterturnedon.1075.5.42Dcloud-shockinteractionInthissection,wetestthe2Dcloud-shockinteractionproblem,whichwetestedinSec-tion4.5.5.Theinitialconditionsweregivenin(4.79)(ˆ;ux;uy;uz;p;Bx;By;Bz)=8>>>>>>>><>>>>>>>>:(3:86859;11:2536;0;0;167:345;0;2:1826182;2:1826182)ifx<0:05;(10;0;0;0;1;0;0:56418958;0:56418958)ifx>0:05,r<0:15,(1;0;0;0;1;0;0:56418958;0:56418958)otherwise,wherer=q(x0:25)2+(y0:5)2.Theinitialmagneticpotentialweusewasgivenin(4.79),Az=8>>><>>>:2:1826182x+0:080921431ifx0:05;0:56418958xifx0:05:InFigures5.12and5.13,weshowresultsobtainedbyrunningthesimulationona256256uniformmeshthatcoversthephysicaldomain(x;y)2[0;1][0;1].WenotethattheHLLDgivesbetterresolutionofshocksandothercomplexfeatures.AsimilarphenomenonisalsoobservedinFigures5.15and5.16,whereweshowresultsobtainedbyrunningthesimulationona256256non-uniformmeshdeterminedbyx=(32˘)cos(ˇ+(12)ˇ=4)+3cos(ˇ=4);y=(32˘)sin(ˇ+(12)ˇ=4)+0:5;(5.93)where(˘;)2[0;1][0;1].108ln(density).kBk.Pressure.Figure5.12:2Dcloud-shockinteraction.AlternativeWENOwithLax-Friedrichs256256uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon.ln(density).kBk.Pressure.Figure5.13:2Dcloud-shockinteraction.AlternativeWENOwithHLLD256256uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon.109Figure5.14:Sketchofthenon-uniformmeshusedincloud-shockproblem.Actualmeshis256256.A3232meshisshownhereforclarity.110ln(density).kBk.Pressure.ln(density).kBk.Pressure.Figure5.15:2Dcloud-shockinteraction.AlternativeWENOwithLax-Friedrichs256256non-uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon.Thesecondrowshowstheregionnearthecenterofthecomputationaldomain.111ln(density).kBk.Pressure.ln(density).kBk.Pressure.Figure5.16:2Dcloud-shockinteraction.AlternativeWENOwithHLLD256256non-uniformmesh.Schlierenplotsatt=0:06.Constrainedtransportandpositivity-preservinglimiterturnedon.Thesecondrowshowstheregionnearthecenterofthecomputationaldomain.112Chapter6ConclusionsandfutureworkInthisdissertationweproposedtwohighorderdischemesfortheidealmag-netohydrodynamics.Theschemeissingle-stagesingle-step.ThebaseschemeusesaWENOmethodwithaLax-WtimediscretizationthatisbasedonthePicardintegralformulationofhyperbolicconservationlaws.Adiscretedivergence-freeconditiononthemagneticisenforcebyusinganunstaggeredconstrainedtransportmethodthatevolvesavectorpotentialalongsidetheconservedquantitiesonthesamemeshastheconservedvariables.ThisvectorpotentialisevolvedwithamoversionofaLax-WWENOmethodthatwasoriginallydevelopedforHamilton-Jacobiequations.Thisallowsustonon-oscillatoryderivativesforthemagneticTofurtherenhancetherobustnessofourscheme,alimiterisaddedtopreservethepositivityofthedensityandpressure.UnlikeourprevioussolversthatarebasedonSSP-RKtimestepping[23,22],thissolverdoesnotrequiretheuseofintermediatestages.Thisreducesthetotalstoragerequiredforthemethod,andmayleadtomoretimplementationforanAMRsetting.Moreover,weneedonlyapplyoneWENOreconstructionpertimestepforthevariables,whereasthethirdandfourth-ordersolversin[23]requirethreeandten,respectively,WENOreconstructionspertimestep.Numericalresultsshowthatourschemehastheexpectedhigh-orderaccuracyforsmoothproblems,andiscapableofsolvingsomeverystringenttest113problems.ThesecondschemeisbasedonanalternativeformulationoftheWENOscheme.Inthisformulation,thenumericalisTaylorexpandedinspace.Theleading(lowestorder)termintheexpansionisapproximatedbyapplyingaRiemannsolveronone-sidedapproximationstotheconservedquantitiesonhalfgridpoints.Theseone-sidedapproxi-mationsareobtainedbyWENOinterpolation.ThehigherordertermsintheexpansionisapproximatedbyacombinationofcentralerencesandanadditionallimiterthatisbasedonthesmoothnessindicatorsforWENOinterpolation.Similartoourscheme,anunstaggeredconstrainedtransportmethodisusedtocontrolthedivergenceerrorofthemagneticThemagneticpotentialisevolvedusingaWENOschemeoriginallydesignedfortheHamilton-Jacobiequationsinthesamewaywedidin[23,22].Alsosimilartoourscheme,apositivity-preservinglimiterisappliedtofurtherenhancetherobustnessofthescheme.Thisschemeisalsoappliedoncurvilinearmeshes,wherewechangetheidealMHDequationstoanewhyperbolicsystemofconservationlawsunderthecurvilinearcoordinatesandchangethemagneticpotentialevolutionequationtoanewequationunderthecurvilinearcoordinates.NumericalresultsshowsthatwhenweusetheHLLDRiemannsolverinthebasescheme,theresolutionofshocksandothercomplexfeaturesisbetterthanwhatwewouldgetbyusingtheLax-FriedrichswhichiscomparabletotheresultswegotbyusingtheclassicalWENOschemewithLax-Friedrichssplitting.Thenumericalresultsalsodemonstratedtherobustnessofourscheme,onbothuniformandnon-uniformmeshes.Thisdissertationservesasthestartingpointforthefollowingdirectionsofpossiblefurtherexplorations.1.Single-stagesingle-stepschemesbasedonthealternativeformulationoftheWENO114scheme.Theworkin[51]suggestedthatthealternativeformulationcanbemadeintoaLax-Wendrtypesingle-stagesingle-stepschemebyapplyingtheCauchy-KovalevskayaproceduresimilartowhatwedescribedinChapter4.SuchaschemewouldhaveanevensmallerestencilthanwhatwehaveinChapter4.2.Treatmentofmorecomplexphysicalboundaryconditions.Theperiodicandzerothorderextrapolationweusedinournumericalexamplesareusefulinareassuchassimulationsofgalaxialformationsandlarge-scalestructuresofmattersintheuniverse.However,inareassuchassolarsystemphysicsandthermonuclearfusion,morecomplexboundaryconditionsareneededanditisanon-trivialtasktoimplementtheminahighorderescheme.3.Incorporatethesingle-stagesingle-stepschemeintoAMRframeworks.Thisisanimportantmotivationforworkinthisdissertation,aswellastheworkin[23,22].115BIBLIOGRAPHY116BIBLIOGRAPHY[1]D.S.Balsara,Second-order-accurateschemesformagnetohydrodynamicswithdivergence-freereconstruction,Astrophys.J.Suppl.S.,151(2004),p.149.[2],Divergence-freereconstructionofmagneticandWENOschemesformagne-tohydrodynamics,JournalofComputationalPhysics,228(2009),pp.5040{5056.[3]D.S.Balsara,Self-adjusting,positivitypreservinghighorderschemesforhydrody-namicsandmagnetohydrodynamics,J.Comput.Phys.,231(2012),pp.7504{7517.[4]D.S.Balsara,C.Meyer,M.Dumbser,H.Du,andZ.Xu,imple-mentationofADERschemesforEulerandmagnetohydrodynamicalonstructuredmeshes{speedcomparisonswithRunge-Kuttamethods,J.Comput.Phys.,235(2013),pp.934{969.[5]D.S.Balsara,T.Rumpf,M.Dumbser,andC.-D.Munz,highac-curacyADER-WENOschemesforhydrodynamicsanddivergence-freemagnetohydrody-namics,J.Comput.Phys.,228(2009),pp.2480{2516.[6]D.S.BalsaraandC.-W.Shu,Monotonicitypreservingweightedessentiallynon-oscillatoryschemeswithincreasinglyhighorderofaccuracy,J.Comput.Phys.,160(2000),pp.405{452.[7]D.S.BalsaraandD.Spicer,Maintainingpressurepositivityinmagnetohydrody-namicsimulations,J.Comput.Phys.,148(1999),pp.133{148.[8]D.S.BalsaraandD.S.Spicer,AstaggeredmeshalgorithmusinghighorderGo-dunovtoensuresolenoidalmagneticinmagnetohydrodynamicsimulations,J.Comput.Phys.,149(1999),pp.270{292.[9]J.W.BanksandW.D.Henshaw,Upwindschemesforthewaveequationinsecond-orderform,J.Comput.Phys.,231(2012),pp.5854{5889.[10]T.J.Barth,Numericalmethodsforgasdynamicsystemsonunstructuredmeshes,inAnintroductiontorecentdevelopmentsintheoryandnumericsforconservationlaws,Springer,1999,pp.195{285.[11]J.A.Bittencourt,FundamentalsofPlasmaPhysics,SpringerNewYork,NewYork,NY,2004.[12]J.P.BorisandD.L.Book,Flux-correctedtransport.I.SHASTA,atransportalgorithmthatworks[J.Comput.Phys.11(1973),no.1,38{69],J.Comput.Phys.,117135(1997),pp.170{186.WithanintroductionbyStevenT.Zalesak,Commemorationofthe30thanniversaryofJ.Comput.Phys.[13]W.Boscheri,M.Dumbser,andD.S.Balsara,High-orderADER-WENOALEschemesonunstructuredtriangularmeshes|applicationofseveralnodesolverstohydrodynamicsandmagnetohydrodynamics,Internat.J.Numer.MethodsFluids,76(2014),pp.737{778.[14]J.U.BrackbillandD.C.Barnes,ThectofnonzerorBonthenumericalsolutionofthemagnetohydrodynamicequations,J.Comput.Phys.,35(1980),pp.426{430.[15]M.BrioandC.C.Wu,Characteristicfortheequationsofmagnetohydrody-namics,inNonstrictlyhyperbolicconservationlaws(Anaheim,Calif.,1985),vol.60ofContemp.Math.,Amer.Math.Soc.,Providence,RI,1987,pp.19{23.[16]L.Chacon,Anon-staggered,conservative,schemefor3Dimplicitex-tendedmagnetohydrodynamicsincurvilineargeometries,Comput.Phys.Commun.,163(2004),pp.143{171.[17]Y.Chen,G.Toth,andT.I.Gombosi,Aderenceschemeforhyperbolicequationsonblock-adaptivecurvilineargrids,J.Comput.Phys.,305(2016),pp.604{621.[18]Y.Cheng,F.Li,J.Qiu,andL.Xu,Positivity-preservingDGandcentralDGmethodsforidealMHDequations,J.Comput.Phys.,238(2013),pp.255{280.[19]A.J.Christlieb,X.Feng,D.C.Seal,andQ.Tang,Ahigh-orderpositivity-preservingsingle-stagesingle-stepmethodfortheidealmagnetohydrodynamicequations,J.Comput.Phys.,316(2016),pp.218{242.[20]A.J.Christlieb,Y.Guc˘lu,andD.C.Seal,ThePicardintegralformulationofweightedessentiallynonoscillatoryschemes,SIAMJ.Numer.Anal.,53(2015),pp.1833{1856.[21]A.J.Christlieb,Y.Liu,Q.Tang,andZ.Xu,Highorderparametrizedmaximum-principle-preservingandpositivity-preservingWENOschemesonunstructuredmeshes,J.Comput.Phys.,281(2015),pp.334{351.[22],Positivity-preservingenceweightedENOschemeswithconstrainedtransportforidealmagnetohydrodynamicequations,SIAMJ.Sci.Comput.,37(2015),pp.A1825{A1845.118[23]A.J.Christlieb,J.A.Rossmanith,andQ.Tang,Finiteenceweightedessentiallynon-oscillatoryschemeswithconstrainedtransportforidealmagnetohydro-dynamics,J.Comput.Phys.,268(2014),pp.302{325.[24]P.Colella,Multidimensionalupwindmethodsforhyperbolicconservationlaws,J.Comput.Phys.,87(1990),pp.171{200.[25]W.DaiandP.R.Woodward,Onthedivergence-freeconditionandconservationlawsinnumericalsimulationsforsupersonicmagnetohydrodynamical,Astrophys.J.,494(1998),p.317.[26]W.DaiandP.R.Woodward,Asimpleenceschemeformultidimen-sionalmagnetohydrodynamicalequations,J.Comput.Phys.,142(1998),pp.331{369.[27]A.Dedner,F.Kemm,D.Kroner,C.-D.Munz,T.Schnitzer,andM.We-senberg,HyperbolicdivergencecleaningfortheMHDequations,J.Comput.Phys.,175(2002),pp.645{673.[28]W.D.D'haeseleer,W.N.G.Hitchon,J.D.Callen,andJ.L.Shohet,ToroidalFluxCoordinates,SpringerBerlinHeidelberg,Berlin,Heidelberg,1991,pp.116{155.[29]A.Dubey,K.Antypas,A.C.Calder,C.Daley,B.Fryxell,J.B.Gal-lagher,D.Q.Lamb,D.Lee,K.Olson,L.B.Reid,P.Rich,P.M.Ricker,K.M.Riley,R.Rosner,A.Siegel,N.T.Taylor,K.Weide,F.X.Timmes,N.Vladimirova,andJ.ZuHone,EvolutionofFLASH,amulti-physicssimulationcodeforhigh-performancecomputing,Int.J.HighPerform.C.,28(2014),pp.225{237.[30]C.R.EvansandJ.F.Hawley,Simulationofmagnetohydrodynamic{aconstrainedtransportmethod,Astrophys.J.,332(1988),pp.659{677.[31]S.A.FalleandS.S.Komissarov,Ontheinadmissibilityofnon-evolutionaryshocks,J.PlasmaPhys.,65(2001),pp.29{58.[32]H.FengandJ.Wang,Observationsofa2!3typeinterplanetaryintermediateshock,SolarPhys.,247(2008),pp.195{201.[33]M.FeyandM.Torrilhon,Aconstrainedtransportupwindschemefordivergence-freeadvection,inHyperbolicproblems:theory,numerics,applications,Springer,Berlin,2003,pp.529{538.[34]T.A.GardinerandJ.M.Stone,AnunsplitGodunovmethodforidealMHDviaconstrainedtransportinthreedimensions,J.Comput.Phys.,227(2008),pp.4123{4141.119[35]S.K.Godunov,Aencemethodfornumericalcalculationofdiscontinuoussolu-tionsoftheequationsofhydrodynamics,Mat.Sb.(N.S.),47(89)(1959),pp.271{306.[36]Godunov,S.K.,Zabrodin,A.V.,Ivanov,M.Ya.,Kraiko,A.N.,andProkopov,G.P.,Numericalsolutionofmultidimensionalproblemsofgasdynamics,Izdat.\Nauka",Moscow,1976.[37]S.V.Golovin,NaturalcurvilinearcoordinatesforidealMHDequations.non-stationarywithconstanttotalpressure,PhysicsLettersA,375(2011),pp.283{290.[38]T.I.Gombosi,K.G.Powell,andD.L.DeZeeuw,Axisymmetricmodelingofcometarymassloadingonanadaptivelyrdgrid:Mhdresults,J.Geophys.Res.{Space,99(1994),pp.21525{21539.[39]S.Gottlieb,C.-W.Shu,andE.Tadmor,Strongstability-preservinghigh-ordertimediscretizationmethods,SIAMRev.,43(2001),pp.89{112(electronic).[40]K.F.Gurski,AnHLLC-typeapproximateRiemannsolverforidealmagnetohydrody-namics,SIAMJ.Sci.Comput.,25(2004),pp.2165{2187(electronic).[41]A.Harten,B.Engquist,S.Osher,andS.R.Chakravarthy,Uniformlyhigh-orderaccurateessentiallynonoscillatoryschemes.III,J.Comput.Phys.,71(1987),pp.231{303.[42]A.Harten,P.D.Lax,andB.vanLeer,OnupstreamencingandGodunov-typeschemesforhyperbolicconservationlaws,SIAMRev.,25(1983),pp.35{61.[43]C.Helzel,J.A.Rossmanith,andB.Taetz,Anunstaggeredconstrainedtransportmethodforthe3Didealmagnetohydrodynamicequations,J.Comput.Phys.,230(2011),pp.3803{3829.[44],Ahigh-orderunstaggeredconstrained-transportmethodforthethree-dimensionalidealmagnetohydrodynamicequationsbasedonthemethodoflines,SIAMJ.Sci.Com-put.,35(2013),pp.A623{A651.[45]W.D.Henshaw,Ahigh-orderaccurateparallelsolverforMaxwell'sequationsonoverlappinggrids,SIAMJournalonScienComputing,28(2006),pp.1730{1765.[46]P.Janhunen,ApositiveconservativemethodformagnetohydrodynamicsbasedonHLLandRoemethods,J.Comput.Phys.,160(2000),pp.649{661.[47]S.C.Jardin,Reviewofimplicitmethodsforthemagnetohydrodynamicdescriptionofmagneticallycdplasmas,J.Comput.Phys.,231(2012),pp.822{838.120[48]A.JeffreyandT.Taniuti,Non-linearwavepropagation,MathematicsinScienceandEngineering,NewYork:AcademicPress,1964,1(1964).[49]G.-S.JiangandC.-W.Shu,implementationofweightedENOschemes,J.Comput.Phys.,126(1996),pp.202{228.[50]G.-S.JiangandC.-c.Wu,Ahigh-orderWENOenceschemefortheequationsofidealmagnetohydrodynamics,J.Comput.Phys.,150(1999),pp.561{594.[51]Y.Jiang,C.-W.Shu,andM.Zhang,AnalternativeformulationofenceweightedENOschemeswithLax-Wendrtimediscretizationforconservationlaws,SIAMJ.Sci.Comput.,35(2013),pp.A1137{A1160.[52],Free-streampreservingenceschemesoncurvilinearmeshes,MethodsAppl.Anal.,21(2014),pp.1{30.[53]M.G.KivelsonandC.T.Russell,IntroductiontoSpacePhysics,CambridgeUniversityPress,1995.[54]P.LaxandB.Wendroff,Systemsofconservationlaws,Comm.PureAppl.Math.,13(1960),pp.217{237.[55]R.J.LeVeque,NumericalMethodsforConservationLaws,LecturesinMathematicsETHZurich,auserVerlag,Basel,seconded.,1992.[56],FiniteVolumeMethodsforHyperbolicProblems,CambridgeTextsinAppliedMathematics,CambridgeUniversityPress,Cambridge,2002.[57]F.Li,L.Xu,andS.Yakovlev,CentraldiscontinuousGalerkinmethodsforidealMHDequationswiththeexactlydivergence-freemagnetic,J.Comput.Phys.,230(2011),pp.4828{4847.[58]S.Li,AnHLLCRiemannsolverformagneto-hydrodynamics,J.Comput.Phys.,203(2005),pp.344{357.[59]S.LiandH.Li,Amoderncodeforsolvingmagnetohydrodynamicsorhydrodynamicequations,tech.rep.,LosAlamosNationalLab.,2003.[60]X.-D.Liu,S.Osher,andT.Chan,Weightedessentiallynon-oscillatoryschemes,J.Comput.Phys.,115(1994),pp.200{212.[61]P.LondrilloandL.D.Zanna,High-orderupwindschemesformultidimensionalmagnetohydrodynamics,TheAstrophysicalJournal,530(2000),p.508.121[62]P.McCorquodaleandP.Colella,Ahigh-ordermethodforcon-servationlawsonlocallyrdgrids,Commun.Appl.Math.Comput.Sci.,6(2011),pp.1{25.[63]A.MignoneandP.Tzeferacos,Asecond-orderunsplitGodunovschemeforcell-centeredMHD:theCTU-GLMscheme,J.Comput.Phys.,229(2010),pp.2117{2138.[64]T.MiyoshiandK.Kusano,Amulti-stateHLLapproximateRiemannsolverforidealmagnetohydrodynamics,J.Comput.Phys.,208(2005),pp.315{344.[65]S.OsherandC.-W.Shu,High-orderessentiallynonoscillatoryschemesforHamilton-Jacobiequations,SIAMJ.Numer.Anal.,28(1991),pp.907{922.[66]J.Qiu,WENOschemeswithLax-WendrtypetimediscretizationsforHamilton-Jacobiequations,J.Comput.Appl.Math.,200(2007),pp.591{605.[67]J.Qiu,M.Dumbser,andC.-W.Shu,ThediscontinuousGalerkinmethodwithLax-Wendrtypetimediscretizations,Comput.MethodsAppl.Mech.Engrg.,194(2005),pp.4528{4543.[68]J.QiuandC.-W.Shu,FiniteenceWENOschemeswithLax-Wendretimediscretizations,SIAMJ.Sci.Comput.,24(2003),pp.2185{2198.[69]P.L.Roe,ApproximateRiemannsolvers,parametervectors,andenceschemes,J.Comput.Phys.,43(1981),pp.357{372.[70]M.RosenbluthandM.Bussac,MHDstabilityofspheromak,Nucl.Fusion,19(1979),p.489.[71]J.A.Rossmanith,Anunstaggered,high-resolutionconstrainedtransportmethodformagnetohydrodynamic,SIAMJ.Sci.Comput.,28(2006),pp.1766{1797.[72]D.D.Schnack,LecturesinMagnetohydrodynamics,vol.780ofLectureNotesinPhysics,Springer-Verlag,Berlin,2009.[73]D.C.Seal,Q.Tang,Z.Xu,andA.J.Christlieb,Anexplicithigh-ordersingle-stagesingle-steppositivity-preservingenceWENOmethodforthecompress-ibleEulerequations,J.Sci.Comput.,(2015),pp.1{20.[74]C.Shen,J.-M.Qiu,andA.Christlieb,AdaptivemeshrbasedonhighorderenceWENOschemeformulti-scalesimulations,J.Comput.Phys.,230(2011),pp.3780{3802.[75]Y.Shen,G.Zha,andM.A.Huerta,E-CUSPschemefortheequationsofidealmagnetohydrodynamicswithhighorderWENOscheme,J.Comput.Phys.,231(2012),pp.6233{6247.122[76]C.-W.Shu,Essentiallynon-oscillatoryandweightedessentiallynon-oscillatoryschemesforhyperbolicconservationlaws,inAdvancednumericalapproximationofnonlinearhyperbolicequations(Cetraro,1997),vol.1697ofLectureNotesinMath.,Springer,Berlin,1998,pp.325{432.[77],Highorderweightedessentiallynonoscillatoryschemesforconvectiondominatedproblems,SIAMRev.,51(2009),pp.82{126.[78]C.-W.ShuandS.Osher,implementationofessentiallynon-oscillatoryshock-capturingschemes,J.Comput.Phys.,77(1988),pp.439{471.[79],implementationofessentiallynon-oscillatoryshock-capturingschemes,II,J.Comput.Phys.,83(1989),pp.32{78.[80]A.Stegmeir,D.Coster,O.Maj,andK.Lackner,Numericalmethodsfor3dtokamaksimulationsusingaeindependentgrid,Contrib.Plasm.Phys.,54(2014),pp.549{554.[81]D.P.Stern,Theartofmappingthemagnetosphere,J.Geophys.Res.-Space,99(1994),pp.17169{17198.[82]J.M.Stone,T.A.Gardiner,P.Teuben,J.F.Hawley,andJ.B.Simon,Athena:anewcodeforastrophysicalMHD,Astrophys.J.Suppl.S.,178(2008),p.137.[83]K.TakahashiandS.Yamada,Regularandnon-regularsolutionsoftheRiemannprobleminidealmagnetohydrodynamics,J.PlasmaPhys.,79(2013),pp.335{356.[84],ExactRiemannsolverforidealmagnetohydrodynamicsthatcanhandlealltypesofintermediateshocksandwaves,J.PlasmaPhys.,80(2014),pp.255{287.[85]V.A.TitarevandE.F.Toro,ADER:arbitraryhighorderGodunovapproach,J.Sci.Comput.,17(2002),pp.609{618.[86]E.F.Toro,Riemannsolversandnumericalmethodsfordynamics,Springer-Verlag,Berlin,thirded.,2009.Apracticalintroduction.[87]E.F.Toro,R.C.Millington,andL.A.M.Nejad,TowardsveryhighorderGodunovschemes,inGodunovmethods(Oxford,1999),Kluwer/Plenum,NewYork,2001,pp.907{940.[88]E.F.Toro,M.Spruce,andW.Speares,RestorationofthecontactsurfaceintheHLL-Riemannsolver,ShockWaves,4(1994),pp.25{34.[89]M.Torrilhon,UniquenessconditionsforRiemannproblemsofidealmagnetohydro-dynamics,J.PlasmaPhys.,69(2003),pp.253{276.123[90]G.Toth,TherB=0constraintinshock-capturingmagnetohydrodynamicscodes,J.Comput.Phys.,161(2000),pp.605{652.[91]K.Waagan,ApositiveMUSCL-Hancockschemeforidealmagnetohydrodynamics,J.Comput.Phys.,228(2009),pp.8609{8626.[92]C.Wang,X.Dong,andC.-W.Shu,ParalleladaptivemeshrmethodbasedonWENOenceschemeforthesimulationofmulti-dimensionaldetonation,J.Comput.Phys.,298(2015),pp.161{175.[93]Z.Xu,Parametrizedmaximumprinciplepreservinglimitersforhighorderschemessolvinghyperbolicconservationlaws:one-dimensionalscalarproblem,Math.Comp.,83(2014),pp.2213{2238.[94]A.L.Zachary,A.Malagoli,andP.Colella,Ahigher-orderGodunovmethodformultidimensionalidealmagnetohydrodynamics,SIAMJ.Sci.Comput.,15(1994),pp.263{284.[95]S.T.Zalesak,Fullymultidimensionalorrectedtransportalgorithmsfor,J.Comput.Phys.,31(1979),pp.335{362.[96]X.ZhangandC.-W.Shu,Onpositivity-preservinghighorderdiscontinuousGalerkinschemesforcompressibleEulerequationsonrectangularmeshes,J.Comput.Phys.,229(2010),pp.8918{8934.[97]U.Ziegler,Acentral-constrainedtransportschemeforidealmagnetohydrodynamics,J.Comput.Phys.,196(2004),pp.393{416.124