Self-Gravitational Magnetohydrodynamics with Adaptive Mesh Refinement for Protostellar Collapse

23 Pages·2007·9.12 MB·English

PASJ:Publ.Astron.Soc.Japan59,905–927,2007October25
(cid:1)c 2007.AstronomicalSocietyofJapan.
Self-Gravitational Magnetohydrodynamics with Adaptive Mesh Refinement
for Protostellar Collapse
TomoakiMATSUMOTO
FacultyofHumanityandEnvironment,HoseiUniversity,Fujimi,Chiyoda-ku,Tokyo102-8160
[email protected]
(Received2006September5;accepted2007June25)
Abstract
D
A new numerical code, called SFUMATO, for solving self-gravitational magnetohydrodynamics (MHD) o
w
problems using adaptive mesh refinement (AMR) is presented. A block-structured grid is adopted as the grid of n
lo
the AMR hierarchy. The total variation diminishing (TVD) cell-centered scheme is adopted as the MHD solver, a
d
e
withhyperboliccleaningofthedivergenceerrorofthemagneticfieldalsoimplemented. TheMHDsolverexhibits d
asecond-orderaccuracyinconvergencetestsoflinearizedMHDwaves.Theself-gravityissolvedusingamultigrid fro
m
method composed of (1) a full multigrid (FMG)-cycle on the AMR hierarchical grids, (2) a V-cycle on these h
grids, and (3) an FMG-cycle on the base grid. The multigrid method exhibits spatial second-order accuracy, fast ttp
s
convergence,andscalability. ThenumericalfluxesareconservedbyusingarefluxingprocedureinboththeMHD ://a
c
solverandthemultigridmethod. Severaltestsareperformedandtheresultsindicatethatthesolutionsareconsistent a
d
withpreviouslypublishedresults. em
Keywords:hydrodynamics—magnetohydrodynamics:MHD—methods: numerical—ISM:clouds—stars: ic
.o
formation up
.c
o
m
/p
1. Introduction the latter approach, all blocks contain the same number of a
s
cells, but the physical sizes of the blocks are different. This j/a
Protostellarcollapseisoneofthemostimportantprocesses approachhasalsobeenadoptedbymanycodes: e.g.,FLASH rtic
le
in star formation; it exhibits a high dynamic range in density (Fryxelletal.2000;Banerjee&Pudritz2006),andNIRVANA -a
b
and spatial dimensions. Adaptive mesh refinement (AMR) (Ziegler 2005). This approach offers advantages including s
is a powerfultechnique forperformingnumericalsimulations relativelysimplealgorithmsforrefinement,parallelization,and trac
requiring a high dynamic range in the mesh schemes. Local vectorization. Thispaperalsofollowsthislatterapproach. In t/5
9
high-resolution is realized by employing grids of differing particular, thisapproachenablestheAMRcodetoimplement /5/9
resolutions. The finer grids are inserted and their location a full multigrid (FMG)-cycle in the multigrid method for 05
changed according to given refinement criteria. Following self-gravity. This scheme is an extension of the multigrid /14
2
the introduction of the fundamental principles of AMR by method for the nested grid approach (Matsumoto & Hanawa 8
8
Berger and Oliger (1984) and Berger and Colella (1989), 2003a)totheAMRgrid. 47
b
this technique is becoming widely used in astrophysical There are several ways to treat MHD, particularly in the y
simulations. treatment of the r (cid:2)B term, where B denotes magnetic gu
e
In star formation, along with self-gravity, the magnetic field. These approaches may be categorized as (1) the st o
field also plays an important role. Both self-gravity and constrained transport method with a staggered grid, (2) n
0
magnetohydrodynamics(MHD) are therefore implemented in the projection scheme with a Poisson solver, and (3) the 9
A
recentAMRcodes(seeareviewbyKleinetal.2007). eight-wave formulation (see To´th 2000; Balsara & Kim 2004 p
Existing self-gravitational MHD AMR codes are generally foracomparisonoftheseschemes). ril 2
0
basedonCartesiangrids;althoughblock-structuredgridstend Recently, Dedneretal.(2002) proposed an alternative 19
to be most often frequently, alternative approaches also exist, scheme for cleaning the r (cid:2)B term. This scheme has
e.g., RAMSES (Fromang et al. 2006). In block-structured a formulation that is similar to, but distinct from, the
AMR,numericalcellsarerefinedintheunitoftheblock,and eight-waveformulation;twoadditionalwavestransferther(cid:2)B
a block is itself an ordinal uniform grid. Block-structured error isotropically at a given speed, and the waves decay at
grids are divided into two categories: patch-oriented grids agivenrate. Asaresult,theerrorispropagatedindependently
and self-similar blocks. In the former approach, a block has of the gas motion, and so is diluted. By contrast, using the
variable number of cells, and the shape of the block is also eight-waveformulation, asingleadditionalwavetransfersthe
changedinthecourseofrefinement. Alargeblockcontaining error at the gas velocity. The error is always propagated
many cells and a small block containing just a few cells can downstream of the gas motion, and can become stagnated at
co-exist. This approach originated in the study ofBergerand the shock wave and the center of the collapsing cloud. The
Colella(1989),andhasbeenadoptedinmanycodes: e.g.,the schemeofDedneretal.(2002)isthereforeadoptedhere.
AMR code Orion of Berkeley (Truelove et al. 1998), Enzo The author has previously developed a self-gravitational
(Norman & Bryan 1999), and RIEMANN (Balsara 2001). In MHD code with a nested grid (Matsumoto & Hanawa 2003a,
906 T.Matsumoto [Vol.59,
2003b; Matsumoto et al. 2004; Machida et al. 2004). In the state, the component corresponding to E in equation (1) and
nestedgridapproach,thegridsarerefinedadaptively,whilethe vectors(3)–(5)(theeighthcomponent)areexcluded,andP is
sizes and positions of the sub-grids are fixed. The problems expressedasafunctionof(cid:2).
to which this method may be applied have therefore been
restricted to those problems whose objects that require fine 3. Discretization
grid resolution are located in the center of the computational
domain,e.g.,Machidaetal.(2006). Equations (1) and (2) are solved by a difference scheme
A self-gravitational MHD code with AMR has now, based on the finite-volume approach. The computational
however, been developed, and this method is presented in domain is divided into cells, each of size ∆x(cid:4)∆y (cid:4)∆z.
this paper. In sections 2–4, the governing equations to be Each cell is labeled by .i;j;k/, the indices of the cell
solved,therulesofdiscretization,andtherefinementalgorithm in the x, y, and z-directions, respectively. The location
are presented. In section 5, the present methods for treating of the cell center is indicated by the position vector r .
i;j;k
the hydrodynamics and MHD of the problem are described, The conservative variables U, the source term S, and the D
o
while in section 6, the methods of the multigrid iteration for gravitational potential Φ are defined at the cell center, i.e., wn
solving self-gravity are described. In section 8, the results of Ui;j;k W=U.ri;j;k/,Si;j;k W=S.ri;j;k/,andΦi;j;k W=Φ.ri;j;k/. loa
d
severalnumericaltestsare presented. Thepaperisconcluded The numerical fluxes F , F , and F are defined at the e
x y z d
insection9. cell surfaces with normals in the x, y, and z-directions, fro
resp(cid:3)ectively. For conv(cid:4)enience, the notation of Fx;i˙1=2;j;k W= m h
2. GoverningEquations Fx ri;j;k˙.∆x=2/xO isintroduced,wherexO denotestheunit ttp
s
vectorinthex-direction. Weintroducethefollowingnotation ://a
Theequationsofself-gravitational,idealhydrodynamicsand todescribethespatialdifferences, c
a
d
MHD can be expressed in a conservative form with a source Qi+1;j;k(cid:3)Qi;j;k em
term, @xQi+1=2;j;k = ∆x ; (6) ic.o
@U +@Fx +@Fy +@Fz =S; (1) @ Q = Qi+1=2;j;k(cid:3)Qi(cid:3)1=2;j;k; (7) up.c
@t @x @y @z x i;j;k ∆x om
togret2hΦer=wi4th(cid:1)GPo(cid:2)i;sson’sequation (2) The dif@fe2xrQeni;cje;ks=inQthi+e1;yj;-k(cid:3)an2d∆Qzxi-;jd2;ikre+ctQioin(cid:3)s1;ja;rke: expressed(8in) /pasj/artic
le
where U is a vector of conservative variables, Fx, Fy, and asimilarmanner. -a
b
F are numerical fluxes, and S is the source term vector. In s
eqzuation(2),Φ,G,and(cid:2)denotethegravitationalpotential,the 4. GridRefinement trac
gravitationalconstant,andthedensity,respectively. t/5
9
ForidealMHDwithself-gravity,thevectorsinequation(1) A self-similarblock-structured gridis adopted. Eachblock /5
areexpressedby consistsofN (cid:4)N (cid:4)N cells,whereN ,N ,andN denote /90
x y z x y z 5
(cid:1) (cid:2) the number of cellsin the x, y, and z-directions respectively, /1
U = (cid:2);(cid:2)vx;(cid:2)vy;(cid:2)vz;Bx;By;Bz;(cid:2)E T; (3) withallblockshavingthesamenumberofcells. Thenumber 428
0 (cid:2)v 1 ofcellsinsideablockisfixedincourseofcalculation,butthe 847
Fx =BBBBBBBBB@ .(cid:2)E+(cid:2)vPx2++(cid:2)(cid:2)PjvvB(cid:3)xxv+jvvvx2yzzB=jBB(cid:3)(cid:3)y8x(cid:1)j(cid:3)BB02+/x=xxvvBB8vyx(cid:1)Byz(cid:3)B==x(cid:3)4Bz4(cid:1)(cid:1)B.x2B=4(cid:2)(cid:1)v/=4(cid:1) CCCCCCCCCA; (4) wcNctgheerhyuillIdllisf=cwiihsnsseiNosdtlsimhfztdh-heeose=siwdemtci8hnfceieflelieialnlnsrprl.sstafihsrdlaAigeiestenuipfistrsiefecsgnybhu1ddleor.aiiemvnc.Tkigradehteoiifiescdnnnaeditunlhmisamteooegbngr8raetreimrcfidcho-nrloiifeetlfeddvcreeaibillonl.lbsnotloT,coiskhct8sehks,ee-cbsathltnboraildultcdockNectscuvkxerealerr=ilndyes by guest on 09 April 20
x x (see figure 1). The cell width of the child cells is half of the 1
(cid:1) (cid:2) 9
S = 0;(cid:2)g ;(cid:2)g ;(cid:2)g ;0;0;0;(cid:2)g(cid:2)v T; (5) parent cell width, that is the refinement ratio is fixed at two
x y z
in thisimplementation (insome AMRshavingpatch-oriented
wherev=.v ;v ;v /T representsvelocity,B=.B ;B ;B /T gridsanyrefinementratioofthepowerof2canbeused). The
x y z x y z
represents the magnetic field, g = .g ;g ;g /T = (cid:3)rΦ coarsestgrid-levelislabeled`=0(thebasegrid),andthefinest
x y z
represents gravity, E = jvj2=2+.(cid:3) (cid:3)1/(cid:3)1P=(cid:2)+jBj2=8(cid:1)(cid:2) grid-level is labeled ` = `max. The `-th grid-level has a 2`
isthetotalenergy,andP representspressure. ThevectorsF timeshigherspatialresolutionthanthecoarsestgrid-level.The
y
andF areobtainedbyrotatingthecomponentsinF bythe block-structured grid is managed by an octree structure; the
z x
right-handrule. parent (coarse) block is linked with eight fine (child) blocks.
For ideal hydrodynamics with self-gravity, the governing Moreover,ablockislinkedwithitsneighboringblocks.These
equationsareobtainedbysettingB=0andomittingthe5thto linklistsarereconstructedeverytimethegridsarerefined.
the7thcomponentsofequation(1)andvectors (3)–(5). In the construction of a child block, the conservative
Barotropic and isothermal equations of state are also variables of the parent cells UH are interpolated to obtain
I;J;K
implemented in the AMR code. In these equations of those of the child cells Uh . For this purpose linear
i;j;k
No.5] Self-GravitationalMHD-AMR 907
5. HydrodynamicsandMHD
5.1. BasicSolvers
The governing MHD equations (1), which include as
aspecialcasethehydrodynamicsequations,aresolvedonthe
block-structuredgriddescribedinsection4. Methodsusedon
ordinal uniform grids can be applied to the block-structured
grid if boundary conditions are properly specified for
each block. A monotone upstream-centered scheme for
conservationlaws(MUSCL)approachandpredictor-corrector
method are adopted here for integration with respect to time
in order to achieve a second-order accuracy in space and
time (e.g., Hirsch 1990), and an unsplit approach rather than D
o
afractionaltimestepapproachisadopted. w
n
Fig.1. Schematicdiagramofthegridrefinementprocess. Thethick ThenumericalfluxisobtainedusingthelinearizedRiemann load
linesrepresenttheboundariesoftheblocksandthethinlinesrepresent solver; the solvers are based on the schemes of Roe (1981) ed
theboundariesofthecells.Eachblockconsistsof83cellsinthisfigure. and Fukuda and Hanawa (1999). For MHD, the hyperbolic fro
Thehatchedblockisrefinedinto8childblocks. divergence cleaning of Dedner et al. (2002) is adopted m
for reducing r (cid:2)B. According to Dedner et al. (2002), http
s
equation(cid:1)s(3)–(5)arethenmodifiedasfollow(cid:2)s: ://a
c
interpolationwithaslopelimiterisadopted, U = (cid:2);(cid:2)vx;(cid:2)vy;(cid:2)vz;Bx;By;Bz;(cid:2)E; T; (12) ade
m
Uhi;j;k=UHI;J;K+rUH (cid:2).rhi;j;k(cid:3)rHI;J;K/; (9) 0 (cid:2)vx 1 ic.o
wthheecreenrtehi;rjs;koafndthreHIc;Jh;Kilddeannodtepathreenptocseitlilos,nrveescpteocrtsivienldyi.caTtinhge BBBB (cid:2)vx2+(cid:2)Pvx+vyjB(cid:3)jB2=xB8(cid:1)y=(cid:3)4(cid:1)Bx2=4(cid:1) CCCC up.com
gradientinsidetheparentcellisslope-limitedaccordingto BB (cid:2)vxvz(cid:3)BxBz=4(cid:1) CC /pa
0 (cid:5) (cid:6) 1 Fx =BB CC; sj/a
B minmod(cid:5)@xUHI+1=2;J;K;@xUHI(cid:3)1=2;J;K(cid:6) C BB vxBy(cid:3)vyBx CC rtic
B C B (cid:3)v B +v B C le
rUH =B@ minmod(cid:5)@yUHI;J+1=2;K;@yUHI;J(cid:3)1=2;K(cid:6) CA; @ .(cid:2)E+P +jBj2z=8x(cid:1)/v x(cid:3)Bz .B(cid:2)v/=4(cid:1) A -ab
x x s
minmod @ UH ;@ UH c2B tra
z I;J;K+1=2 z I;J;K(cid:3)1=2 h x c
(10) (13) t/59
0 1 /5
where minmod.(cid:2);(cid:2)/ denotes the minmod function. Note that B (cid:2)g (cid:3).r(cid:2)0B/B =4(cid:1) C /905
this interpolation conserves the conservative variables in the B x x C /1
refiZnementprocedure, Z BBB (cid:2)(cid:2)ggy(cid:3)(cid:3)..rr(cid:2)(cid:2)BB//BBy==44(cid:1)(cid:1) CCC 4288
B z z C 4
UH.r/dr = Uh.r/dr; (11) S =BB 0 CC; (14) 7 b
ΩHI;J;K ΩHI;J;K BB 0 CC y gu
B 0 C e
whereΩHI;J;K denotesthezoneofaparentcellwhosecenteris @ (cid:2).g(cid:2)v/(cid:3)B(cid:2).r /=4(cid:1) A st o
locatedatrHI;J;K. (cid:3).ch=cp/2 n 09
The refinement algorithm is based on that of Berger and A
Cfoollloewllain(g1p9r8o9c)e,dwurheesretogcroidnsstroufctlegvreilds`oafrleevreelfi`n+ed1u:sing the wchhiesreth e wisavaesscpaelaerd,paontedntcipalisprtohpeadgaamtinpginagdraivteerogfenthceewerarvoer., pril 2
0
Note that this modification requires additional components in 1
9
1. A block of grid-level` is markedif the cells inside the the basic equation, and a total of nine waves are solved for.
blocksatisfyarefinementcriterion. This approach is similar to the eight-wave formulation. In
2. Ifblocksofgrid-level`+2exist,thenthecorresponding the method of Dedner et al. (2002), two additional waves
blocksof`arealsomarked. transfer the r (cid:2)B error isotropically at a speed of c , and
h
3. Blocks adjacent to marked blocks are also marked, so the waves decay at a rate of c . Similarly, in the eight-wave
p
that the grids of level `+1 properly nest the grids of formulation,anadditionalwavetransfersthedivergenceerror
level`+2. atagasvelocityofv.
4. Newblocksoflevel`+1areconstructedaschildblocks In the predictor step, Un is updated to Un+1=2 by a half
ofthemarkedblocks. i;j;k i;j;k
timestep,i.e.,
5. Blocksof`+1areremovediftheirparentblocksarenot (cid:5)
∆t
marked. Un+1=2=Un (cid:3) @ Fn +@ Fn
i;j;k i;j;k 2 x x;i(cid:6);j;k y y;i;j;k
These procedures are called in ascending order of grid-level,
from`max(cid:3)1to0. +@zFnz;i;j;k(cid:3)Sni;;j(cid:5);k ; (15)
908 T.Matsumoto [Vol.59,
where the superscript n denotes the time level, and The source term in the corrector step is estimated at the time
∆t =tn+1(cid:3)tn. The numerical flux Fn, which is defined at levelt =tn+1=2,
thecellboundaries,iscalculatedusingtheprimitivevariables,
(cid:1) (cid:2) Sn+1=2=
T i0;j;k 1
Q= (cid:2);v ;v ;v ;B ;B ;B ;P; : (16)
x y z x y z 0
The numerical flux has second-order spatial accuracy due to BB (cid:2)n+1=2gn+1=2(cid:3).r(cid:2)B/n+1=2Bn+1=2=4(cid:1) CC
B i;j;k x;i;j;k i;j;k x;i;j;k C
theMUSCLextrapolationoftheamplitudesoftheeigenmode B (cid:2)n+1=2gn+1=2(cid:3).r(cid:2)B/n+1=2Bn+1=2=4(cid:1) C
oftheMHDwaves,LU,whereLdenotesthematrixoftheleft B i;j;k y;i;j;k i;j;k y;i;j;k C
eigenvector (see the appendix in Fukuda & Hanawa 1999 for BBB (cid:2)in;+j;k1=2gzn;+i;j1;=k2(cid:3).r(cid:2)B/ni;+j;k1=2Bzn;+i;j1;=k2=4(cid:1) CCC;
details)1. B 0 C
B C
ThesourcetermSn;(cid:5) iscalculatedusing B 0 C
0 i;j;k 1 BB 0 CC
BB (cid:2)n gn(cid:3)1=2(cid:3).r(cid:2)0B/n Bn =4(cid:1) CC @ (cid:2)in;+j;k1=2.gni;+j;k1=2(cid:2)vni;+j;k1=2/(cid:3)0Bni;+j;k1=2(cid:2).r /ni;+j;k1=2=4(cid:1) A Down
BBB (cid:2)iin;;jj;;kkgxyn;;(cid:3)ii;;1jj=;;kk2(cid:3).r(cid:2)B/ini;;jj;;kkBxyn;;ii;;jj;;kk=4(cid:1) CCC (20) loade
Sni;;j(cid:5);k =BBBBB (cid:2)in;j;kgzn;(cid:3)i;1j=;k2(cid:3).r(cid:2)0B/ni;j;kBzn;i;j;k=4(cid:1) CCCCC; gwnh+er1e=2(cid:2)isin;o+j;bk1=ta2inaenddbyvnis;+jo;lk1v=2ingarPeoiosbsotani’nseedqufarotimon(U2)in;+ja;tk1=th2,ehaanldf d from h
BBB@ (cid:2)in;j;k.gni;(cid:3)j;1k=2(cid:2)vni;j;k/(cid:3)00Bni;j;k(cid:2).r /ni;j;k=4(cid:1) CCCA timi;j(cid:1)e;rks2teΦp(cid:2)an+s1f=o2ll=ow4s(cid:1),G(cid:2)n+1=2: (21) ttps://aca
0 i;j;k i;j;k de
m
(17) Poisson’sequationissolvedbythemultigridmethod,asshown ic
in section 6, and the multigrid method prepares not only the .ou
where .r (cid:2)B/i;j;k = @xBx;i;j;k +@yBy;i;j;k +@zBz;i;j;k, and cell-centeredΦ,butalsothecell-surfacedg. Gravityatthecell p.c
o
the magnetic field defined at the cell surfaces, Bx;i˙1=2;j;k, is centergn+1=2 isobtainedbyaveragingthevaluesofgravityat m
.orbt ai/nie;jd;kf=rom.@xth ei;nj;ikn;t@hyc oim;j;pk;o@nze nti;jo;kf/FTxi;si˙o1b=t2a;ji;nke.dSfirmomilartlhye, thecellsiu;jr;kfaces0, 1 /pasj/a
esfiFovfztlahu;liut;jitao;oktne˙sd1ies=vs2uee.psneTathdrhaaectesonlmfyionpbltlohyonwoceospn,metsrpaotoofnreFsnpxtl;iio˙tfti1nt=h2g;e,j;iksn,ouwFrhcyie;ci;hjte˙trh1m=e2;(fk1o,4rma)naidsl gni;+j;k1=2= 21B@ gggynxnzn;+;+;i+ii;+;j11j1==+;1=k22=21+2=;12j=;;kk2+++gggxyznnn;;;+++iii;;(cid:3)jj111;===1(cid:3)k=222(cid:3)12=1;j2=;;2kk CA: (22) rticle-abstra
n+1=2= n+1=2;(cid:5)exp"(cid:3)∆2t (cid:7)ccph(cid:8)2#; (18) tNthhoeetreecfetolhlr-aectecgnotieanrtceidadecΦselw[lsisetuehrfetahqceueactieisonontrbsatl(a4idn2ief)f–de(r4be4yn)c]te.heoEfdqiΦuffaetariseonnfca(er2ao2sf) ct/59/5/905
where n+1=2;(cid:5) is the ninth component of Un+1=2 solved by theinsideofablockisconcerned. /142
Similar to the predictor step, the ninth component of the 8
equation(15). Thefreeparametersc andc arerelatedtothe 8
h p sourcetermisevaluatedbyoperatorsplitting, 47
timNeo-mteatrhcahtinggra,vaintydgarne(cid:3)d1=e2sclargibsebdyinhaslufbaseticmtieonst5e.p2..Thisslight " (cid:7)c (cid:8)2# by g
i;j;k n+1= n+1;(cid:5)exp (cid:3)∆t h ; (23) ue
laggingisexpectedtohaveanegligibleeffectontheaccuracy c s
inthe predictor step, andthegravityinthe previouscorrector p t on
stepcanbereusedinthispredictorstep(Trueloveetal.1998). where n+1;(cid:5) is the ninth component of Un+1 obtained by 09
Thisavoidsanadditionalcallofthemultigridmethodtosolve equation(19). Ap
Poisson’s equation, significantly reducing the computational 5.2. TimeMarching ril 2
costsofthismethod. 01
9
In the corrector step, a spatially second-order numerical The time-marching represented by equations (15) and (19)
fluxFn+1=2 isobtainedbyapplyingMUSCL extrapolationto proceedsinunitsatthegrid-level. Thecodeisequippedwith
the amplitudes of the eigenmodes, which are converted from two modesof time-marching: an adaptiveand a synchronous
Un+1=2. Using this flux, Un is updated to Un+1 by a full time-stepmode. Whichmodeisuseddependsonwhetherthe
i;j;k i;j;k
timestep, gasisnon-self-gravitationalorself-gravitational. Theadaptive
(cid:5) time-stepmodeisappropriatefornon-self-gravitationalgases,
Uni;+j;k1=Uni;j;k(cid:3)∆t @xFnx+;i;1j=;k2+@yFny+;i;1j=;k2 and is based on the method of Berger and Colella (1989). In
(cid:6)
this mode, a coarser grid has a longer time step than a finer
+@zFnz;+i;j1=;k2(cid:3)Sni;+j;k1=2 : (19) grid. In contrast, the synchronous mode is appropriate for
self-gravitationalgases,andeverygrid-levelhasthesametime
step. This is because evolution on the fine grid affects the
1 Forthecaseofthehydrodynamics,theMUSCLextrapolationisapplied detached coarse grid immediately due to self-gravity, and so
totheprimitivevariables, Qinthepredictorandcorrectorstepsinthis thesametimestepmusttobechosenforeverygrid-level(see
implementation.
No.5] Self-GravitationalMHD-AMR 909
In both cases, c is constant across all the grid-levels, and
h
satisfiestheCFLconditionateverygrid-level.
Thedampingratec isobtainedfromc ,
p h
c2 =0:18Lc ; (30)
p h
where L is a scale length of a problem, and the coefficient
of 0.18 is chosen according to Dedner et al. (2002). This
coefficient is a free parameter specifying the ratio of the
damping rate and propagation speed of r(cid:2)B. We confirmed
that this selection worked well in many calculations (e.g.,
Machida et al. 2005a, 2005b) and in the numerical tests
presentedinthispaper.
D
o
5.3. BoundaryConditionsforGhostCells w
n
Fig.2. Adaptive time-step scheme. The thick arrows denote time lo
Each block has ghost cells overlapping with the adjacent a
sotfetpism,ew-mhialercthhienga.ssociatednumbers,(cid:1)1 ,(cid:1)2 ,...,(cid:1)9 ,indicatetheorder blocks. Figure 3 shows the interface between the coarse and ded
fineblocksinthex-direction. Twoghostcellsarepreparedin fro
m
each direction, as indicated by the gray circles in the figure, h
discussion of Truelove et al. 1998). Note that the adaptive sincetheMUSCLextrapolationrequirestwocellsonboththe ttp
s
time-stepmodecouldalsobeusedforaself-gravitationalgasat leftandrightsidesofacellboundarywherenumericalfluxis ://a
theexpenseofthefirst-ordertemporalaccuracyofthescheme. obtained. ca
d
Figure 2 shows the order in which the grid-levels proceed A boundary condition is imposed on the ghost cells each e
m
for the adaptive time-step mode schematically. The numbers time before the predictor and corrector steps proceed. When ic
.o
associated with the thick arrows denote the order of a block is adjacent to another block of the same grid-level, u
p
time-marching. Coarser grid-levels precede finer grid-levels. data is exchanged between the two blocks by a simple copy. .co
Thefinegrid-levelundergoesseveralsub-cyclesuntilthetime Whenablockisadjacenttoacoarseblock,dataisinterpolated m/p
level of the fine grid-level is synchronized with that of the spatially. In addition, data is interpolated temporally if an as
coarse grid-level. The time step of a finer grid-level, ∆th, is adaptive time-step mode is used. The procedure for the j/a
givenby boundaryconditionsisasfollows: rticle
∆th=∆tH2(cid:3)n; (24) 1. Adjacentcoarseblocksareidentified. -ab
n=min˚m2Nj∆tH2(cid:3)m(cid:6)∆tChFL(cid:9): (25) 2. Cthoengsheorvsattcievlelsvoafriaafiblneesbolonctkhaerecoinatresrepocelalltsedotveemrlpaoprpailnlyg stract/5
wCFhLereco∆ndtChitFioLndaetntohteesfintheegtriimd-elesvteelp,acnadlc∆ultaHteddedniroetcetslythbeytitmhee tporoccoeinducirdeeiswoitmhitthteedtiimfethleevseylnocfhtrhoenfiounsetgirmide--lsetveepl.mTohdies 9/5/90
5
stepatthecoarsergrid-level. Notethat∆th isfixedat∆tH=r isused. /1
4
in Berger and Colella (1989), where r denotes a refinement 3. Thetemporallyinterpolatedvariablesareconvertedinto 2
8
ratio.Ontheotherhand,inthemethodpresentedhere∆thcan primitive variables. They are interpolated spatially to 84
7
beequalto∆tH=2nforn=0;1;2;(cid:2)(cid:2)(cid:2),if∆thsatisfiestheCFL coincidewiththegridpointsoftheghostcellsofthefine b
y
condition. Inthesynchronoustime-stepmode,acommontime grid-level. Theinterpolatedvariablesarecopiedintothe g
u
step∆t isusedatallgrid-levels,andisgivenby ghostcells. es
˚ (cid:9) 4. Adjacentblocksinthesamegrid-levelareidentified. t o
∆t =0(cid:6)m`(cid:6)i`nmax ∆tC`FL ; (26) 5. Tsihmeplvyarcioapbileeds ionftoctehllesgohvoesrtlcaeplplsin.g the ghost cells are n 09 A
wcohnedrietio∆ntCa`tFgLridde-lneovteels`.the time step calculated by the CFL In the temporal interpolation, the quadratic interpolation is pril 20
The solver described in subsection 5.1 includes two performed on UH;n, UH;n+1=2, and UH;n+1. Note that UH;n 19
parameters, c and c , and these are related to the mode of and UH;n+1 have of second-order temporal accuracy while
h p
thetime-marching.Thewavespeedc isobtainedas UH;n+1=2 is of first-order temporal accuracy. Using these
h
variables,UH.t/intn(cid:6)t (cid:6)tn+1isinterpolatedas,
h
(cid:7) (cid:8)
ch=CFL∆t`=0; (27) (cid:1) (cid:2) t(cid:3)tn 2
UH.t/= UH;n+1(cid:3)U.1/;H;n+1
where CFL denotes the CFL number, and ∆t`=0 denotes the tn+1(cid:3)tn
timestepatgrid-level`=0. Fortheadaptivetime-stepmode, .tn+1(cid:3)t/UH;n+.t(cid:3)tn/U.1/;H;n+1
˚ (cid:9) + ; (31)
h=min ∆x`;∆y`;∆z` ; (28) tn+1(cid:3)tn
`=0 where U.1/;H;n+1 denotes the conservative variables of
whileforthesynchronoustime-stepmode, first-orderaccuracyatthetimeleveltn+1,definedas,
˚ (cid:9)
h= min ∆x`;∆y`;∆z` : (29) U.1/;H;n+1=2UH;n+1=2(cid:3)UH;n: (32)
`=`max
910 T.Matsumoto [Vol.59,
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
Fig.3. Theinterfacebetweenacoarseblockandafineblock.Thegridpointsoftheghostcellsforthefineblockaredenotedbygraycircles,whilethe ic
realgridpointsaredenotedbyfilledcircles.Thearrowsdenotethenumericalfluxesthroughtheinterfacebetweenthefineandcoarseblocks. .o
u
p
.c
o
m
/p
P a
iinnstAeidrpeostlhlaoetpioceen-ll.limI.nIitfi;eJgd;uKreg/r3ai,sdtiheevenagtlruaaidsteiednautdsooifnptgtheedprfiomritivtheevarsipaabtlieasl ainccjOfio;kOgr;dnuirFneghx4;t;nio,(cid:3)1t=h2e;jO;krOe∆flSuxh∆ingth;pn.rocMedourreeeuxppdliactietsly,UaHs;(cid:5)intdoicaUteHd sj/article-a
rQH =0BB@ minmod(cid:5)(cid:5)@yQ@xHIQ;J+HI(cid:3)1=1=22;K;J;;K@yQHI;J(cid:3)1=2;K(cid:6)(cid:6) 1CCA: UH =UH;(cid:5)(cid:3)∆V1H Xn suXrfaceFhx;n∆t!h;n∆Sh bstract/59/5
/9
minmod @ QH ;@ QH 0
z I;J;K+1=2 z I;J;K(cid:3)1=2 (33) (cid:3)FHx ∆tH∆SH ; (34) 5/142
8
8
where 4
It may be noted that equation (33) describes a simple X 7 b
interpolation in the normal direction, with a slope-limited ∆th;n=∆tH; (35) y g
u
interpolationinthetransversedirections. Xn es
5.4. RefluxingattheInterfaces ∆Sh=∆SH; (36) t on
0
surface 9
RefluxingwasintroducedbyBergerandColella(1989),and A
has often been used in standard implementations of AMRs. and∆VH denotesthevolumeofthecoarsecell. pril 2
This technique maintains consistency between the fine and 0
1
coarse grid-levels during time-marching, and ensures that 6. MultigridMethodforSelf-Gravity 9
conservationlaws,e.g.,massconservation,aresatisfied. 6.1. MultigridCycles
Asshowninfigure3,UH isupdatedusingFH ,
I;J;K x;I+1=2;J;K The multigrid method is widely used in many AMRs for
while Uh is updated using Fh;n through several
i;jO;kO x;i(cid:3)1=2;jO;kO solving Poisson’s equation. In many AMRs, a solution
sub-cycles until the fine grid-level catches up with the converges by means of the V-cycle on the grids of AMR
coarse grid-level, where jO 2 Œj;j + 1(cid:4), kO 2 Œk;k + 1(cid:4), hierarchy, and the FMG-cycle on a uniform base grid. By
and n denotes the index of the sub-cycles of the fine contrast, the multigrid method presented here uses not only
grid-level. Obviously, FH ∆SH∆tH should equal V-cycles, but also FMG-cycles on the grids of the AMR
P x;I+1=2;J;K
jO;kO;nFhx;;ni(cid:3)1=2;jO;kO∆Sh∆th;n to ensure that the conservation himiepralermcheyn.tedAninFMtheG-acpypclreoaocnh tohfe Mhieartsaurcmhoictaol agnriddsHiasnaalwsoa
laws are satisfied, where ∆SH = ∆yH∆zH and ∆Sh = (2003a),andthesamestrategyisadoptedhere.
∆yh∆zh. In refluxing, UHI;J;K is re-calculated, taking Figure5showsthegridsusedinthemultigridcyclesforthe
accountofthedifferencebetweenFH ∆SH∆tH and caseN =N =N =8. Forconvenience,agridislabeledby
x;I+1=2;J;K x y z
No.5] Self-GravitationalMHD-AMR 911
6.2. Smoothing
The red-black Gauss-Seidel iteration is adopted as
a smoothing procedure. Conventionalred-blackGauss-Seidel
iterationcanbe appliedonlyto a uniform grid. Thisiteration
scheme is therefore modified so that it can be applied to
composite grids, in which grids with different resolutions
co-exist.
The discretization of Poisson’s equation on the composite
gridsisnowdescribed.Poisson’sequationcanbeexpressedas
asetoftwoequations,
Fig.4. Refluxingprocedureattheinterfacebetweenfineandcoarse
cells. r(cid:2)g=(cid:3)4(cid:1)G(cid:2); (39)
g=(cid:3)rΦ: (40) Do
w
n
Theseequationsarediscretizedas lo
a
whΩerm`e`awnidthmd`e2noŒt0e;t`hmeagxr(cid:4)id-laenvdelomfth2eŒ0A;MmRmahxi(cid:4)e;rarchy(a3n7d) @gxgx;i;j;k+=@y(cid:3)gy@;i;Φj;k+@zgz;;i;j;k =(cid:3)4(cid:1)G(cid:2)i;j;k; ((4412)) ded from
the coarsening level, respectively. For example, the hatched x;i+1=2;j;k x i+1=2;j;k h
grid in figure 5e is labeled as Ω2, where the grid of `=1 is gy;i;j+1=2;k =(cid:3)@yΦi;j+1=2;k; (43) ttps
schooawrsneniendfibgyuraefa5catoisrm22.=T0h,ewcoh1ailresethnaintgofletvheelgorfidasllsthhoewgnridins gz;i;j;k+1=2=(cid:3)@zΦi;j;k+1=2; (44) ://acad
figure5bism=1,irrespectiveofthegrid-level`. Inaddition where gx;i+1=2;j;k, gy;i;j+1=2;k, and gz;i;j;k+1=2 are the em
acompositegridisdefinedas componentsofgravitydefinedonthecellsurfaces. ic.o
[ [ Considering thefine cellsadjacentto theinterfacebetween u
ΩOmW=Ωm0 (cid:2)(cid:2)(cid:2) Ωm`max: (38) thefineandcoarsecellsinfigure7,thepotentialontheghost p.com
Forexample,thecompositegridsshowninfigures5a,5b,and cell,ΦB isrequiredinordertoobtaingxh;i(cid:3)1=2;j;k. Thepotential /pa
5careexpressedasΩO0,ΩO1,andΩO2,respectively. ontheghostcellisgivenby, sj/a
An overview of the cycles used by the present numerical 10Φh +8Φ(cid:5)(cid:3)3Φh rtic
O ΦB = i;j;k i+1;j;k; (45) le
method is now presented. First, the solution on Ωm for 15 -a
b
0fiFIngM(cid:6)uoGrtehm-secry5(cid:6)cwal–emo5Orcimdssaa,xgntidhvceeo6nnagv.rbeiydrTghiemsOesmcmuoanaaxdrxs=eiemrnluFeomdMg2uGcmno-taciilnrysct(cid:1)ehlNneeisxn,n;guNamslyeb;vsNeherlozwo(cid:2)infn(cid:3)tthhi1nee. Φ(cid:5)= 9ΦHI;J;K+3(cid:1)ΦHI;J;K(cid:3)1+16ΦHI;J(cid:3)1;K(cid:2)+ΦHI;J(cid:3)1;K(cid:3)1;(46) stract/59/5/9
cellsperblockisdecreaseduptotwoinatleastonedirection. 05
At the bottom of the FMG-cycle, which is marked by V in where ΦB is obtained by quadratic interpolation in the /14
figure 6a, the solution then converges on Ωm for 0 (cid:6) ` (cid:6) x-direction,andΦ(cid:5) isobtainedbybilinearinterpolationinthe 28
fi`mguaxresan5dd–m5f)=. TmOypmiacxalluyn,doenrlyVo-cnyeclietesra(ttihoenhoa`ftcthheedV-gcryicdlseiins yto–tzheplianntee.rfTacheesqautiasdfireasticthientneercpeoslsaatiroyncionntdhietiodnirsecfotirosnon-coarlmleadl 847 by
sufficient,although,forreference,twocyclesoftheV-cycleare conservativeinterpolation(seeTrottenbergetal.2001). Using gu
e
illustrated in figure 6b. At the bottom of the V-cycle, which this procedure, the components of gravity of the fine cells st o
is marked by B in figure 6b, the solution on the coarsened adjacenttocoarsecellsareobtained. n
0
basegridΩm0 form(cid:7)mOmax convergesunderFMG-cycles, as Forthe coarsecellsadjacentto fine cells, thegravityatthe 9 A
shown in figures 5g–5i and 6c. Finally, at the bottom of the interfaceisgivenbysummingthegravityonthecorresponding p
FMG-cycleonthebasegrid,whichismarkedbyEinfigure6c, cellsurfacesofthefinecells. Fortheexampleoffigure7,such ril 2
0
anexactsolutionisgivenaccordingtotheboundaryconditions, acoarsegravityisgivenby, 19
sincethereisonlyonecellinthecomputationaldomain. jX+1Xk+1
O 1
NotethatthecomputationatΩ0 oftheFMG-cycleaccounts gH = gh (47)
x;I+1=2;J;K 4 x;i(cid:3)1=2;jO;kO
for most of the computational tiOme of the multigrid method, jO=jkO=k
because the number of cells in Ω0 dominates all others. The
number of cells decreases by a factor of 1/8 every time the instead of gxH;I+1=2;J;K = @xΦHI+1=2;J;K. This is a similar
grids are coarsened. This indicates that the computation strategy to the refluxing described in subsection 5.4, and
on ΩO0 dominates the overall computational cost, and this ensures that the solution satisfies Gauss’s theorem; when the
normal component of gravity is summed over the surfaces
scheme is scalable to the number of cells in a similar way as
of any cell, the sum equals the mass contained by the cell
the hydrodynamics scheme (see also Matsumoto & Hanawa
multipliedby4(cid:1)G,
2003a).
X
g(cid:2)∆S =(cid:3)4(cid:1)G(cid:2)∆V; (48)
surface
912 T.Matsumoto [Vol.59,
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/p
a
s
j/a
rtic
le
-a
b
s
tra
c
t/5
9
/5
/9
0
5
/1
4
2
8
8
4
7
b
y
g
u
e
s
t o
n
0
9
A
p
ril 2
0
1
9
Fig.5. Schematicdiagramofthecoarseningofgridsinthemultigridmethod.Thecell-boundariesandblock-boundariesaredenotedbythinandthick
lines,respectively. (a–c)CoarseningofgridsintheFMG-cycleontheAMRhierarchy. Thenumberofthecellsperblockdecreasesupto23. (d–f)
V-cycleontheAMRhierarchy. Thesolutionconvergessequentiallyonthehatchedblocks. (g–i)CoarseningofgridsintheFMG-cycleontheuniform
basegrid.
No.5] Self-GravitationalMHD-AMR 913
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
ic
.o
u
p
.c
o
m
/p
a
s
j/a
rtic
le
-a
b
s
tra
c
t/5
9
/5
/9
0
5
/1
4
2
8
8
4
7
b
y
g
u
e
s
t o
n
0
9
A
p
ril 2
0
1
9
Fig.6. Schematicdiagramofthemultigridcycles.Thelinespointingdiagonallydownwardsfromlefttorightdenotetherestrictionoperators,whilethe
linespointingdiagonallyupwardsfromlefttorightdenotetheprolongationoperators. TheSsymbolsdenotethesmoothingoperators. TheVsymbols
inpanel(a)denotetheV-cyclesontheAMRhierarchyshowninpanel(b).TheBsymbolsinpanel(b)denotetheFMG-cyclesontheuniformbasegrid
showninpanel(c).AtthepointdenotedbyE,anexactsolutionisobtainedaccordingtotheboundaryconditions.
914 T.Matsumoto [Vol.59,
D
o
w
n
lo
a
d
e
d
fro
m
h
ttp
s
://a
c
a
d
e
m
Fig.7. Interfacebetweenthefineandcoarsecellsforthemultigridmethod. ic
.o
(cid:1) (cid:2)i up
where∆V denotesthevolumeofthecell,and∆S denotesthe +3 ΦHI;J˙1;K˙1+ΦHI˙1;J;K˙1+ΦHI˙1;J˙1;K ; .com
vectorofthecross-sectionofthecellinthedirectionnormalto (51) /p
thecellsurface.Wealsoconfirmedthatrefluxingthegravityis as
nsuebcseescsatiroynfo8r.5s)e.cond-orderaccuracybyaconvergencetest(see cwehller.ei;.jI;;kJ/;,Kan/dintdhiecasitgesnaofco˙arsdeepceenlldosvoenrlathpepipnagriwtyitohfaifi,nje, j/article
This refluxing of the gravity is adopted only for the andk. -a
b
s
FMG-cycle on the composite grids shown in figure 6a, in The restriction procedure is performed by averaging the tra
which the solution converges simultaneously on Ωm` over the valuesofthefinecells.i;j;k/whichoverlapthecorresponding ct/5
grid-levels `. The refluxing of gravity is not adopted in the coarsecell.I;J;K/, 9
/5
Vgr-icAdy-cclceleov,redlssini,ngcetothethseoludtiisocnretciozantvieorngesofseeqquueanttiioanlly(4o1v)e,r tthhee ΦHI;J;K =IhHΦh= 81Xi+1Xj+1Xk+1ΦhiO;jO;kO: (52) /905/14
Gauss-Seideliterationisexpressedas, iO=ijO=jkO=k 288
4
h2 Theprolongationandrestrictionintroducedaboveisusedin 7 b
Φni;je;wk =Φi;j;k(cid:3) 6 Ri;j;k; (49) commonwithallthemultigridcycles. y gu
e
6.4. FMG-CycleonAMRHierarchy s
Ri;j;k W=@xgx;i;j;k+@ygy;i;j;k+@zgz;i;j;k+4(cid:1)G(cid:2)z;i;j;k t on
= (cid:3)LΦ +4(cid:1)G(cid:2) ; (50) An FMG-cycle based on the standard algorithm of the 09
i;j;k z;i;j;k FMG-cycle for linear equations (see e.g., Press & Teukolsky A
p
LwhdeerneohtedsetnhoeteLsapthlaecciealnlowpiedrtaht,oRr.i;Ej;qkudaetnioonte(s4t9h)eurpedsaidteusalΦ,and 1hi9e9r1a)rchisy,imbpecleamuseenttehdeobnastihceocpoemraptoorssiteargeripdrsepoafretdheoAnMthRe ril 20
i;j;k 1
to Φnew at every iteration. Typically only two iterations are compositegrids. Thesmoothing,prolongation,andrestriction 9
i;j;k
requiredfortheFMG-c(cid:1)ycle,andon(cid:2)eiterationfortheV-cycle. procedures are given by equations (49), (51), and (52),
Hereafter,Φnew=L(cid:3)1 Φ ;(cid:2) referstotheGauss-Seidel respectively. In the smoothing procedure, the refluxing of
i;j;k GS i;j;k i;j;k
the gravity is performed according to equation (47). In
iterationgivenbyequations(49)and(50).
the prolongation, the variables on Ωm are transferred onto
`
6.3. ProlongationandRestriction Ωm(cid:3)1. Similarly, the restriction procedure transfers variables
`
The full-weight prolongation is adopted here, and is given on Ωm onTto Ωm+1. The overlap of coarse grids on fine
` `
by grids, Ωm Ωm , is therefore not taken into account in the
` `+1
FMG-cycle.
Φh = IhΦH
i;j;k H
h 6.5. V-CycleonAMRHierarchy
1
W= 27ΦH +ΦH
64(cid:1) I;J;K I˙1;J˙1;K˙1 (cid:2) In the V-cycle, the multilevel adaptive technique (MLAT)
+9 ΦH +ΦH +ΦH is adopted by using the full approximation scheme (FAS)
I˙1;J;K I;J˙1;K I;J;K˙1