Numerical non-LTE 3D radiative transfer using a multigrid method

14.8 MB·

by  Johan P. Bjørgen

Astronomy&Astrophysicsmanuscriptno.accepted_aa cESO2017
(cid:13)
January9,2017
Numerical non-LTE 3D radiative transfer using a multigrid method
JohanP.Bjørgen1 andJorritLeenaarts1
InstituteforSolarPhysics,DepartmentofAstronomy,StockholmUniversity,AlbaNovaUniversityCentre,SE-10691Stockholm,
Sweden
Received;Accepted
ABSTRACT
Context.3Dnon-LTEradiativetransferproblemsarecomputationallydemanding,andthissetslimitsonthesizeoftheproblemsthat
canbesolved.SofarMultilevelAcceleratedLambdaIteration(MALI)hasbeentothemethodofchoicetoperformhigh-resolution
7
computationsinmultidimensionalproblems.ThedisadvantageofMALIisthatitscomputingtimescalesas (n2),withnthenumber
1 O
ofgridpoints.Whenthegridgetsfiner,thecomputationalcostincreasesquadratically.
0
Aims.Weaimtodevelopa3Dnon-LTEradiativetransfercodethatismoreefficientthanMALI.
2
Methods.Weimplementanon-linearmultigrid,fastapproximationstoragescheme,intotheexistingMulti3Dradiativetransfercode.

n WeverifyourmultigridimplementationbycomparingwithMALIcomputations.Weshowthatmultigridcanbeemployedinrealistic
a problemswithsnapshotsfrom3Dradiative-MHDsimulationsasinputatmospheres.
J Results.Withmultigrid,weobtainafactor3.3-4.5speedupcomparedtoMALI.Withfull-multigridthespeed-upincreasestoafactor

6 6.Thespeedupisexpectedtoincreaseforinputatmosphereswithmoregridpointsandfinergridspacing.
Conclusions.Solving3Dnon-LTEradiativetransferproblemsusingnon-linearmultigridmethodscanbeappliedtorealisticatmo-

] sphereswithasubstantialspeed-up.
R
1. Introduction We therefore need a robust method that can solve the 3D non-
S
LTEproblemmoreefficientlythancurrentlypossible.
.
hMostofourunderstandingofastrophysicalobjectscomesfrom
It is well known that Λ-iteration converges prohibitively
panalysing their spectra. The comparison of models of those
slowly in optically thick non-LTE problems. Accelerated Λ-
-objects to observations therefore requires solving the radiative
o iteration (ALI, Cannon 1973), converges much faster by using
transfer(RT)problemwiththemodelasinput.Inthispaperwe
r anapproximateΛ-operator.For3Donetypicallyusesthediag-
tinvestigateanon-linearmultigridmethodforsolvingthethree-
s onal of Λ-operator as the approximate operator (this choice is
adimensional (3D) non-LTE RT problem in cool stellar atmo-
also known as Jacobi iteration), because it can be easily con-
[spheres,andspecifically,thesolaratmosphere.
structed without the need of computationally expensive matrix
1 One-dimensional non-LTE radiative transfer has been pos- inversions. Other methods exists: Gauss-Seidel (GS) iteration
vsibleformanydecades.Itismorerecentthatmultidimensional (Trujillo Bueno & Fabiani Bendicho 1995) converges twice as
7radiative transfer has become possible. Now that 3D radiation- fast as the diagonal operator, but cannot be parallelized as effi-
0MHD simulations of solar and stellar atmospheres are widely ciently as Jacobi iteration (Tsitsiklis et al. 1988), a serious dis-
6used (Asplund et al. 2003, 2004; Uitenbroek 2006; Trujillo advantage in parallel solvers that are used for 3D applications.
1Bueno & Shchukina 2007; Pereira et al. 2009; Holzreuter & Conjugategradientmethodsareunderdevelopment,buthavenot
0
Solanki 2013; de la Cruz Rodríguez et al. 2013; Pereira et al. yetbeenimplementedinageneral3Dnon-LTERTcode(Pale-
.
12015; Rathore & Carlsson 2015; Amarsi et al. 2016) it has be- tou&Anterrieu2009).
0comeimportanttohaveefficient3Dnon-LTEradiativetransfer
Therefore, multilevel accelerated lambda iteration (MALI)
7methodstocomparethesemodelswithobservations.
(Rybicki & Hummer 1991, 1992) with a diagonal approximate
1
: In the solar chromosphere, full 3D RT is important for operator remains the current standard for 3D non-LTE radia-
vstronglyscatteringlines(e.g.Leenaartsetal.2009,2012,2013a; tive transfer problems for solar applications (Uitenbroek 2001;
XiŠteˇpánetal.2015;Sukhorukov&Leenaarts2016).Thenextso- Leenaarts&Carlsson2009;Šteˇpán&TrujilloBueno2013).
largenerationtelescopes,suchasthe4-meterDKIST,willobtain
r ThedisadvantagewithMALIisthattheconvergenceslows
aadiffractionlimitof 0.03(cid:48)(cid:48) at500nm,correspondingto32km downasthegridspacinggetsfiner(Olsonetal.1986).Asolu-
on the Sun. To be able to compare these high-resolution obser-
tion to this problem was proposed by Steiner (1991), who im-
vations with MHD simulations, 3D radiative transfer is needed
plemented a linear multigrid method for two-level atoms with
evenformuchweakerphotosphericlines(Holzreuter&Solanki
coherentscatteringin1Dand2Dtestatmospheresandfounda
2015; Judge et al. 2015). Some important chromospheric diag-
speed-up of a factor 4 to 40. The idea behind multigrid meth-
nosticlinesareinfluencedbypartialredistributioneffects,such
asLyman-α,Mgiih&k,andCaiiH&K.Theseeffectshavebeen odsisthatJacobiandGauss-Seideliterationquicklysmoothout
thehighspatial-frequencyerrorinthesolution,butdecreasesthe
includeda3DRTcode(Sukhorukov&Leenaarts2016),butat
low-frequencyerroronlyslowly.Afterafewiterationstheerror
the cost of an increase in computational time up to a factor 10.
willthereforebesmooth(i.e.,containonlylowfrequencies).The
problemcanthereforebetransferredtoacoarsergrid,wherethe
Send offprint requests to: J. P. Bjørgen e-mail: lowspatial-frequencyerrorsbecomehigher-frequencyerrors.It-
[email protected] erationsonthecoarsegridquicklydecreasetheseerrors.Acor-
Articlenumber,page1of12
A&Aproofs:manuscriptno.accepted_aa
rectiontothefine-gridsolutioniscomputedonthecoarsegrid, ν0/ν1 ν1+ν2 ν2
k
and interpolated up onto the fine grid, which now has a much
smallerlow-frequencyerror. ν1 ν2 ν1 ν2
This method was expanded by Vath (1994), who applied it k−1
to a 3D problem with a smooth atmosphere. He showed that
the method works, but does not provide a large speed-up on
k−2
small domains in parallel setup. Fabiani Bendicho et al. (1997) ν3 ν3
combined the non-linear multigrid technique with a multilevel ν1+ν2 ν2
k
Gauss-Seideliterationscheme(MUGA).Thesmoothingproper-
tiesofMUGAdonotdeterioratewithincreasinggridresolution, ν ν ν ν
2 2 1 2
whichmakesitan excellentchoicetousewithmultigrid. They k−1
foundaspeed-upofafactor3to32withafive-levelCaiiatom
anda2Datmospherethatwasisothermalintheverticaldirection k−2
andcontainedaweaktemperatureinhomogeneityinthehorizon- νfmg ν3 ν3
tal direction (see Auer et al. 1994, Eq. 18). Léger et al. (2007)
comparedtheperformanceofmultigridagainstGauss-Seidelit-
erationwithsuccessiveoverrelaxationinsmooth2Dprominence Fig.1.SchematicrepresentationofV-cycleandfullmultigridwiththree
grid levels. means interpolating the full approximation, means
models. Šteˇpán & Trujillo Bueno (2013) proved that multigrid restrictingth⇒eapproximationandtheresidual,(cid:100)meansinte→rpolating
works in full 3D non-LTE with a MALI operator and paral-
the correction. The iteration parameters ν and ν are defined in
lel implementation using spatial domain decomposition. They FMG 1,2,3
Section2,theparameterν inSection4.12.
0
used a plane-parallel semi-empirical FAL model C atmosphere
Fontenlaetal.(1993),withasinusoidalhorizontalinhomogene-
ityintemperaturewithanamplitudeof500K. for all frequencies and ray propagation directions. The combi-
All previous studies used relatively smooth model atmo- nation of the equations can be written as a coupled system of
spheres. So far no one has applied non-linear multigrid for non-linearequations:
’reallife’atmospheresproducedbyradiation-MHDsimulations
which are highly inhomogeneous in all atmospheric quantities.
(n)= f. (4)
Thispaperpresentsanon-linearmultigridmethodthatcanhan- W
dlesuchatmospheres. Wenotethatthesearen n equationsthatinprinciple
levels gridpoints
In Section 2 and 3 we describe the method and our com- dependonalln n levelpopulations
levels gridpoints
putationalsetup.InSection4wepresentnumericalconsidera- For a given grid point, if we replace the first rate equation
tions for the 3D radiative transfer with multigrid. In Section 5 withparticleconservation, f lookslike:
we present the speedup we obtain compared to regular MALI
 
n
i2te.rMatieotnh.oSdection6containsourconclusions. f = 0t...ot (5)
0
2.1. Thenon-LTEradiativetransferproblem
IntheMALIschemetheseequationsarepreconditionedand
The time-independent non-LTE radiative transfer problem con-
solvedbyiteration(Rybicki&Hummer1991,1992).
sistsofsolvingthetransportequation
dI
ν =S (n) I , (1) 2.2. Thenon-linearmultigridmethod
dτ (n) ν i − ν
ν i
In this subsection we present how the non-linear multigrid
withI theintensity,τ theopticalpathalongaray,S thesource
ν ν ν method is applied to the radiative transfer problem, i.e. how it
function, and n the atomic level populations, together with the
i isusedtosolveEq.4.
equationsofstatisticalequilibriumforthelevelpopulations:
The non-linear multigrid, fast approximation storage (FAS)
(cid:88) (cid:88) scheme as formulated by Brandt (1977), preserves the mathe-
n P (I ) n P (I )=0, (2)
i ij ν − j ji ν maticalstructureof atthecoarsergrids.Fornon-LTEcodes
W
withP thesumoftheradiativeandcollisionalrates.Wewrote thismeansthatonecanreusetheformalsolver,computationof
ij
the dependencies of τ and S on the level populations and P collisionalrates,etc.,andthattheiterativeMALImethodcanbe
ν ν ij
on the intensity explicitly to emphasize the non-linearity of the usedonthecoarsergridswithonlyminormodifications.
problem.Theradiationfieldcouplesvariousgridpointstogether, Werestatetherateequation:
makingtheproblemalsonon-local.Weassumeinthispaperthat
k(nk)= fk, (6)
thecollisionaltermsandthebackgroundopacityandemissivity W
areconstantandfullydeterminedbytheinputatmosphere. wherekisanindexthatdenotesthegridlevel:k=(cid:96),(cid:96) 1,...,1
Toclosethesystem,wereplaceoneoftherateequationsat where(cid:96)isthenumberofgrids,withk=(cid:96)thefinestgri−dandk=
eachgridpointwithaparticleconservationequation: 1thecoarsestgrid.Weprovideaninitialguessfornandperform
a number (ν ) MALI iterations, which is called pre-smoothing.
n(cid:88)levels 1
n = n. (3) Because MALI iterations act as a smoothing operator, we have
tot i
reducedthehigh-frequencyerrorandobtainedasmoothresidual
i=1
rk:
Thestatisticalequilibriumequationmustbesolvedatallgrid
points,andthetransferequationmustbesolvedatallgridpoints rk = fk k(nk). (7)
−W∗
Articlenumber,page2of12
JohanP.BjørgenandJorritLeenaarts:Numericalnon-LTE3Dradiativetransferusingamultigridmethod
105 105 1016
K] 11001145 3cm]− K] Tne 11001145 3cm]−
e[ y[ e[ 1013 y[
atur 104 1013 nsit atur 1012 nsit
mper 1012 onde mper 104 1011 onde
Te ctr Te 1010 ctr
Temp 1011 Ele 109 Ele
ne
103 108
0.1 2.0
Pre-smoothing
0.0
1.5 CGC
0.1 Post-smoothing
n∞−0.2 n∞ 1.0
/ − /
) )
∞ 0.3 ∞ 0.5
n − n
− 0.4 −
ni − ni 0.0
( 0.5 Pre-smoothing (
−0.6 CGC −0.5
− Post-smoothing
0.7 1.0
− −
-90 11 253 703 1454 2118 2145 2148 2158 0 500 1000 1500 2000 2500 3000 3500 4000
Height[Km] Height[Km]
Fig.2.BehaviouroftherelativeerrorinpopulationatthefinestgridduringvariousstepsinaV-cyclewiththreegrids,fortheFAL-Catmosphere
(left-handpanels)andacolumnextractedfroma3Dradiation-MHDatmosphere(atmosphereModel1,seeSection3).Bothcomputationsare
done for an Hi atom using 10 pre-smoothings, 32 coarse-grid iterations and 25 post-smoothings. The upper panel shows the temperature and
electrondensityintheatmosphere.ThelowerpanelsshowtheerrorfortheHin=2levelafterthepre-smoothing(red),theerrorafterapplication
ofthecoarse-gridcorrection(blue)andtheerrorafterapplicationofthepost-smoothing(purple).Forbothatmospheresthecoarsegridcorrection
hasdecreasedthelowspatial-frequencyerror,butincreasedthehigh-frequencyerror.Thelatterislargelyremovedbythepost-smoothing.
The , denotes that an extra formal solution has been per- griditerations.Therelativecorrection(SeeSection4.6)isthen
∗
formedtogettheradiativeratesconsistentwiththecurrentpop- computedby:
ulationestimate.
Theexactsolutioncanbedescribedbyaddingacorrectionto ek 1 = nakf−te1r−nkb−ef1ore, (11)
tihsethaepnp:roximation:nkexact =nk+ek.Theexactfinegridequation − nbk−ef1ore
wherethesubscriptsbeforeandafterdenotethepopulationsbe-
k(nk+ek)= fk. (8) foreandaftertheν MALIiterations.
W 3
Thecorrectionismappedontotheoperatortothefinergrid
Using the residual equation (Eq. 7) the fine grid equation withaninterpolationoperator k :
(Eq.8)maybewrittenas Ik−1
(cid:16) (cid:17)
nk =nk 1+ k (ek 1) (12)
k(nk+ek)=rk+ k(nk). (9) new old Ik−1 −
W W∗ The interpolation introduces high-frequency errors on the
Allthesequantitiescanbemappedtoacoarsergridk 1usinga finegrid.Byapplyingν2MALIiterationsonthefinegrid(called
restrictionoperator k 1.Wepointoutthatweassum−ethatonly post-smoothing), these high-frequency errors are removed, and
the residual is smooRthk−and not the populations themselves. The wehaveobtainedabetterapproximationtothetruesolutionon
populationsfromthefinegridareusedasaninitialapproxima- the finest grid. To obtain a converged solution, multiple cycles
tionatthecoarsegrid.Thecoarsegridoperator k 1isobtained havetobeperformed.
−
by performing a formal solution at the coarseWgrid. This results The method we just described is called two-grid FAS. The
inthecoarse-grid-equation: sameprocedurecanbeappliedtoEq.10becauseithasthesame
structureasEq.4.Doingthisleadstothegeneralmultigridpro-
cedure. The successive coarsening starting from the finest grid
(cid:16) (cid:17)
k 1(nk 1+ek 1)= k 1(rk)+ k 1 k 1(nk) . (10) and the interpolation steps from the coarsest to the finest grid
W − − − Rk− W∗− Rk− togetherarecalledaV-cycle.
Theright-hand-sidetriestoforcethesolutionofthesystem
tomaintainthefinegridaccuracy.Thecoarsegridequationdoes
2.3. Fullmultigrid
notneedtobesolvedexactly,sinceitislimitedbytheaccuracy
of the right-hand-side. Instead we compute an approximate so- Full multigrid (FMG) can be seen as a method to obtain a bet-
lution by applying a number ν MALI iterations, called coarse ter initial approximation at the finest level. FMG is considered
3
Articlenumber,page3of12
A&Aproofs:manuscriptno.accepted_aa
tobetheoptimalmultigridsolver(Trottenbergetal.2000).We 3.1. Modelatmospheres
explain FMG using a three-grid ((cid:96) = 3) example, as illustrated
Weusethreedifferentmodelatmospheresinthispaper.
inFigure1.InFMG,Eq.4issolvedatthecoarsestgrid:
Thefirstisasnapshotfromtheradiation-MHDsimulationby
Carlssonetal.(2016)computedwiththeBifrostcode(Gudiksen
(cid:16) (cid:17)
k=1 nk=1 = fk=1. (13) etal.2011)takenatt=3850s
W This model extends from the upper convection zone to the
coronainaboxwith504 504 496gridpoints,withahorizontal
The initial population and the right-hand-side can be ob-
× ×
gridspacingof48km,andaverticalgridspacingrangingfrom
tained by restricting it from the finest grid or direct initializa-
19 to 100 km. The model contains strong gradients in temper-
tionatthek=1grid.Afterapplyingaν iterations,aninitial
FMG
ature, velocity and density and has been used before in several
approximationisacquiredatthecoarsestlevel.Thefullapprox-
studies (e.g. de la Cruz Rodríguez et al. 2013; Leenaarts et al.
imationistheninterpolateduptothenextfinergridk=2:
2013a,b; Šteˇpán et al. 2015; Pereira et al. 2015; Loukitcheva
(cid:16) (cid:17)
nk=2 = k=2 nk=1 . (14) etal.2015;Rathore&Carlsson2015;Rathoreetal.2015).We
Ik=1 used this model to see how multigrid behaves in a realistic use
The new approximation is used as the initial population for case.Toobtainanoddnumberofgridpointsintheverticaldirec-
solvingEq.4atthek = 2grid.NowaV-cycleisappliedatthis tion(seeSection4.3),weaddedoneextralayerintheconvection
level, instead of interpolating up to the finest grid k = 3. The zonethroughlinearextrapolationofallquantities,soweobtain
reasonforthisisthattheapproximationatk = 2willhaveboth 504 504 497gridpoints.ThisMHDsnapshotwillbecalled
× ×
low-frequencyandhigh-frequencyerrorsduetotheinterpolation atmosphereModel1.
anddiscretizationerrors.Theseerrorsareefficientlyreducedby WealsouseaBifrostsimulationsnapshotwith768 768
× ×
a V-cycle. When the V-cycle is completed, we interpolate the 768gridpoints(Carlsson&Hansteen,privatecommunication).
population at k = 2 to the finest grid k = 3. Now we have an Inthissnapshotthex-andy-axisareequidistantwithagridspac-
initial approximation at k = 3 that is closer to the true solution ing of 32 km. The vertical axis is non-equidistant, with a grid
thanastraightinitializationatgridk=3. spacing13kminthephotosphereandthechromosphere,andin-
Rememberthatwearenotaimingtogetthefullsolutionfor creasing to 60 km in the convection zone and the corona. The
theradiativetransferproblem,butratherabetterinitialapproxi- physicalextentisthesameasforModel1:24Mm 24Mm
× ×
mationatthefinestgrid.Afterreachingthefinestlevel,theFMG 16.9 Mm. We only used the part of the atmosphere between
isfollowedbyregularFASV-cyclestoobtainthesolution. zmin = 0.5 Mm and zmax = 6 Mm, to cover the formation
height o−f the Caii lines. The grid size is therefore reduced to
InFigure1weshowaschematicrepresentationoftheFAS
768 768 473.Thissimulationhasasimilarmagneticfieldcon-
and full-multigrid-algorithms. Figure 2 shows an illustration of
× ×
figurationasModel1,butwascomputedwithanLTEequation
thebehaviouroftherelativeerrorinasmoothsemi-empiricalat-
ofstateinsteadofthenon-equilibriumequationofstateofModel
mosphereandanon-smoothcolumnfromaradiation-MHDsim-
1.ThisatmospherewillbecalledModel2.Finallywealsouse
ulation, but with otherwise identical setup. Comparing the two
somecomputationsusingtheplane-parallelsemi-empiricalFAL
cases shows that the radiation-MHD simulation column shows
modelCatmosphere(Fontenlaetal.1993),whichwewillrefer
much more high-frequency error after the coarse grid correc-
toasFAL-C.
tion than in FALC. These high frequency errors stem from the
inevitable inability to accurately represent a non-smooth atmo-
sphere on a much coarser grid. We will show in this paper that
3.2. Modelatoms
multigrid can still handle such atmospheres, but at the cost of
anincreaseincomputationaltimeandthusasmallerincreasein We use three different model atoms. Most test computations
computationalspeedcomparedtoMALIthanreportedinearlier were done with a minimalistic three-level Caii atom, contain-
worksforsmoothatmospheres. ing the ground level, the upper level of the Caii K line and the
Caiii ground state. We also use the 6-level Caii atom and the
six-levelhydrogenatomfromCarlsson&Leenaarts(2012).All
3. Computationalsetup bound-boundtransitionsaretreatedincompleteredistribution.
WeimplementedtheFASmultigridmethodincludingFMGinto
theexistingradiativetransfercodeMulti3D(Leenaarts&Carls- 3.3. Convergencecriteriaanderrorestimates
son 2009). This code solves the statistical equilibrium equa-
Tovalidateandcompareourmultigridimplementationwiththe
tionsforoneatomicspecieswithabackgroundcontinuumusing
standardone-gridMALIscheme,wecompareagainstreference
MALI with a diagonal approximate operator (Rybicki & Hum-
solutionscomputedusingtheone-gridMALIscheme.Withthe
mer 1991, 1992). The code used 3D atmospheres on a Carte-
three-level Caii atom we converged the reference solution to
siangridasinput.ItisparallelizedwithMPI1 usingspatialdo-
max δn/n 10 5. For the hydrogen atom we converged to
main decomposition as well as decomposition over frequency. max|δn/n| ≤ 10 −4. The reason for the lower criterium for hy-
Theatomicbound-boundtransitionscanbetreatedincomplete | | ≤ −
drogen was to limit computational expenses: to reach the limit
or partial redistribution (Sukhorukov & Leenaarts 2016). The
forhydrogenweneededroughly300,000CPUhours.
transport equation is solved using short-characteristic (Kunasz
To characterize the absolute error after i iterations we used
&Auer1988)withalinearorHermite-spline(Auer2003;Ibgui
theinfinitynorm:
et al. 2013) interpolation scheme. We use the 24-angle quadra-
(cid:12) (cid:12)
ture (A4 set) from Carlson et al. (1963). All computations pre- (cid:12)n n (cid:12)
sented in this paper are done assuming complete redistribution Ei,∞ =max(cid:12)(cid:12)(cid:12) in− ∞(cid:12)(cid:12)(cid:12), (15)
andwithoutfrequencyparallelization. ∞
wherethen populationsaretakentobethereferencesolutions.
1 MessagePassingInterface Sinceonegr∞idpointcanslowdownorstalltheconvergencewith
Articlenumber,page4of12
JohanP.BjørgenandJorritLeenaarts:Numericalnon-LTE3Dradiativetransferusingamultigridmethod
multigrid,theinfinity-normisthebestchoice.Othernorms,such
asthe2-norm,cangiveafalsepictureoftheconvergenceofthe
multigridmethod.
Themaximumrelativecorrectioninpopulationduringeach
iterationisgivenby
(cid:12) (cid:12)
(cid:12)n n (cid:12)
Ci =max(cid:12)(cid:12)(cid:12) i−n i−1(cid:12)(cid:12)(cid:12). (16)
i
For linear non-LTE radiative transfer problems one can define
thespectralradiusasthelargesteigenvalueoftheamplification
matrixoftheiterationscheme.However,inthenon-linearcase,
this is not possible. Therefore we define the spectral radius op-
erationallyas
E
ρ = i+1. (17)
sr E
i
Itmeasureshowmuchtheerrorisreducedforeachiteration.
For normal use cases there is no reference solution, so the Fig.3.ConvergencebehaviorforMALI(red)andmultigridwiththree
error cannot be computed directly. For MALI iterations it was (blue)andfour(purple)gridsforthethree-levelCaiiatomwithinatmo-
shown (Auer et al. 1994; Fabiani Bendicho et al. 1997) that sphereModel1.Thecomputingtimeisexpressedinunitsofthetime
onecanusethemaximumrelativecorrectiontoestimatetheer- it takes to perform one MALI iteration at the finest grid. The multi-
ror. When the relative correction reaches the asymptotic linear gridcomputationswereperformedwithν1 = 2,ν2 = 2,ν3 = 32,full-
regime, the relation between the error, the spectral radius and weightingrestrictionandtrilinearinterpolation.
therelativecorrectionis:
ρ
E sr C. (18) 4.3. Sub-domainsizesandnumberofgrids
i ≈ 1 ρ i
sr
− Multi3D divides the computational domain into sub-domains.
Becausethereisnoreferencesolution,thespectralradiusises- Theformalsolverrequires5ghostzonesinthehorizontaldirec-
timatedusingρsr Ci+1/Ci.Wetestwhetherthisapproximation tions,andthenumberofinternalpointsinthehorizontaldirec-
≈
isalsovalidformultigriditerationinSec4.13. tionmustbeatleastaslargeasthenumberofghostzones.Inthe
verticaldirectionwehaveonlyoneghostzone.Insingle-gridap-
plications,thesubdomainshavetypicalsizesfrom323to643.In
4. Numericalconsiderations
multigrid,weusethesamedomaindecompositionasinthefinest
Inthefollowingsubsectionswediscussanumberofissuesthat grid,andcoarsenthegridlocaltoeachsub-domain.Thislimits
arerelevantforradiativetransferusingmultigridmethods. thenumberofcoarseningsthatwecando.Fora323subdomain,
wecanachievetwocoarseningto163 and83.A643 subdomain
can be coarsened 3 times. Through various tests we found that
4.1. Gridcoarsening
threegridsgivesusthebestandstableperformance.Usingmore
Weuseasimplecoarseningstrategyforthegridcoordinates,by than three grids can lead to non-convergence or lower perfor-
removingeverysecondgridpointalongthex,yandzaxeswhen mancethanusingthreegrids(seealsoTable1ofSteiner1991).
coarseningfromgridktok 1. InFigure3wedisplaytheresultofsuchatestforthethree-level
− CaiiatomandatmosphereModel1.
4.2. Gridsizes
4.4. Restrictionoperator
Intheverticaldirectionweusefixednon-periodicboundaries.In
the lower boundary, the incoming radiation is set to the Planck We implemented three different restriction operators: injection,
functionatthelocalgastemperature.Intheupperboundary,the half-weightingandfull-weighting(e.g.Trottenbergetal.2000).
incoming radiation is set to zero or to pre-specified values. We Injectionisthemoststraightforwardmethod,asitjustremoves
chosetokeeptheverticalboundariesatfixedheights.Toobtain gridpointswhenmovingaquantityfromafinetoacoarsegrid.
this,werequireanoddnumberofgridpointsintheverticaldi- Half-weighting is a smoothing operation, where a coarse grid
rection at all grid levels. Therefore the number of vertical grid pointvalueisanaverageofthesamepointplusitssixneighbour
pointsN atthefinestlevelmustfulfilthefollowingrelation: gridpointslocatedalongthex,y,andzaxesonthefinegrid.Full-
Z
weighting is similar, but includes all 26 neighbour grid points.
N = A 2(cid:96)+1, (19) Full-weighting has the advantage that it conserves a quantity
z
·
when it is integrated over the entire computational domain. In-
whereAisanintegerand(cid:96)thenumberofgrids.
jection is computationally more efficient than the half and full
Inthehorizontalplaneweuseperiodicboundaryconditions.
weighting methods, but the extra computing time is negligible
BecauseMulti3Drequiresconstantgridspacinginthehorizontal
comparedtothetimespentontheformalsolution.
directionwethusrequirethat N ,thenumberofgridpointsin
x,y WeperformedvarioustestusingboththeFAL-CandModel
thexorydirection,follows
1 atmospheres. We found that injection performs well for the
N = B 2(cid:96), (20) FAL-Catmosphere,butnotforModel1.Figure4showacom-
x,y
· parison between single-grid MALI, three-grid injection, and
withBaninteger. three-grid full-weighting using the three-level Caii atom in
Articlenumber,page5of12
A&Aproofs:manuscriptno.accepted_aa
101 100
MALI1g 32 32 31
× ×
MG3g:fw 62 62 61
× ×
MG3g:inj 96 96 48
100 The×oret×ical
E∞10−1 tFS 10−1
10−2
10−2
0 50 100 150 200 250
1 2 3
Computingtime(MALIiteration)
Grid
Fig.5. Thecomputationalcostforoneformalsolutionateachgridin
Fig. 4. Convergence behavior for MALI (red) and multigrid with in- multigridwithdifferentsub-domainsizes.Wesetthecostofaformal
jection(purple)andfull-weighting(blue)asrestrictionoperators,using solution at the finest grid at 1. The theoretical value is calculated as
thethree-levelCaiiatomandatmosphereModel1.Themultigridcom- (1/8)(cid:96)−k,withkthegridlevel.
putationwasperformedwithν = 2,ν = 2,ν = 32,threegridsand
1 2 3
trilinearinterpolation.
4.7. Theatmosphereatvariousgridlevels
Multigrid requires solving the formal solution at all grid lev-
els. We thus need a method to generate coarse-grid model at-
Model 1. The convergence of multigrid using injection stalls
mospheresfromthefine-gridatmosphere.Wetestedbothinjec-
around E = 0.03, while the full-weighting computation con-
tionandfullweightingtodoso.Ourtestsindicatethatinjection
vergesto∞thereferencesolution.Wethereforechosetousefull-
canleadtothebestconvergencerate,especiallyinsmoothatmo-
weighting for the atomic level populations and residuals for all
sphereslikeFAL-C,butitmightalsocausestallsinconvergence
furthercomputations.
orevendivergence.Full-weightinghasalowerbest-caseconver-
gence, but does not lead to stalls or divergence. The reason for
this difference is that injection more accurately represents the
4.5. Interpolationoperator fine-gridproblem,andsohasbetterbest-caseperformance.But
italsoenhancesgradientsonthecoarsergrid,leadingtostronger
We implemented two interpolation operators: linear and cubic non-linearitiesandlargerprobabilityofstalling.Full-weighting
Hermite(Auer2003;Ibguietal.2013).Thecoefficientsforlin- ontheotherhandleadstocoarse-gridatmospheresthataremore
ear interpolation depend only on the grid spacing and can be differentfromthefine-gridatmosphere,leadingtolowerperfor-
precomputed. A full 3D interpolation uses information of the mance,butalsohassmallergradientsleadingtoincreasedstabil-
8 surrounding grid points. Cubic Hermite interpolation needs a ity.Becausestabilityisusuallythegreatestconcernwhenusing
4 4 4 stencil for a 3D interpolation and requires the com- radiation-MHD atmospheres, we use full-weighting for the at-
× ×
putation of the derivatives of the quantity to be interpolated on mosphereinitializationinourcomputations.
theinner8gridpoints.ThismakescubicHermiteinterpolation
somewhat more involved to implement and more computation-
4.8. Initializingtheradiationfieldatvariousgridlevels
ally expensive. We performed a number of 1D test calculation
using columns from Model 1 and found that, on average, both Multi3D treats background processes using LTE opacities and
interpolationmethodsperformequallywell.Wethereforechose a source function with a thermal part and a coherent-scattering
tousethesimplerandfasterlinearinterpolationforour3Dcom- part.Ateachgridlevel,wethusneedtoperformasmallnumber
putations. ofΛ-iterationstoensurethebackgroundprocessesareinitialized
correctly.BecauseaΛ-iterationtakesasimilartimeasaMALI
iteration we want to minimize the number of Λ-iterations. Em-
4.6. Correction piricallywefoundthatoneΛ-iterationissufficientatthefinest
grid.Atthecoarsergridsweobtainaninitialguessfortheradia-
The radiative transfer problem is highly non-linear, and a level tionfieldbyrestrictingitfromthenextfinergridandperforming
populationcaninprinciplechangebyordersofmagnitudefrom oneortwoΛ-iterations.
one grid point to the next. This poses a problem for the inter-
polationofthecorrectionfromacoarsertoafinergrid(Eqs.11
4.9. Computationalcost
and12)becausetheinterpolationoperationcanintroducealarge
amount of interpolation noise and even lead to the unphysical In this paper we are mainly concerned with the theoretical per-
prediction of negative populations at the fine grid. Both cases formance of multigrid compared to MALI iteration. We there-
destroythegoodconvergencebehaviourofmultigrid.Wefound foreexpresscomputationtimeinunitsofoneMALIiterationat
that interpolating the relative correction (Eq. 11) decreases the thefinestgridinourfigures.Thecomputingtimeofaformalso-
interpolationnoiseandlargelyeliminatestheoccurrenceofneg- lutionatcoarsegridk 1isthen1/8timesthecomputingtime
−
ativepopulationsasaconsequenceofinterpolation. atgridk.
Articlenumber,page6of12
JohanP.BjørgenandJorritLeenaarts:Numericalnon-LTE3Dradiativetransferusingamultigridmethod
However, Multi3D does not follow this theoretical perfect
102
scaling. The reason is that the short-characteristic solver in Zeroradiation
Multi3D requires 5 ghost zones in the horizontal direction and LTEpopulations
101
one ghost zone in the vertical direction. As the number of grid
r
e
points per subdomain gets smaller with increasing coarsening, mb 100
the relative amount of time spent in updating the ghost zones u
n
increases.InFig.5wecomparethetheoreticalmaximumspeed- ng 10−1
upwiththeactualspeedupmeasuredonaHPClusterPlatform hi
ot
3000withIntelXeonE-526002.2GHzcoresusingthree-grids. mo 10−2
At grid k = 2 the performance penalty is a factor 3 for the S
32 32 31subdomain,butinterestinglythe96 96 48sub- 10−3
× × × ×
domainshowsaspeedupbetterthanthetheoreticalvalue.Atthe
k = 1gridtheperformancepenaltyisvisibleforallsubdomain 10−4
sizes that we tested, but with the expected behaviour that the 0 20 40 60 80 100 120 140
largerthesubdomain,thesmallertheperformancepenalty. ν1
4.10. Occurrenceofnegativepopulations Fig.7.DependenceofthesmoothingnumberρL(ν1)(seeEq.22),on
theinitialconditionofthepopulations.Theresultswerecomputedusing
StandardMALIiterationkeepsatomiclevelpopulationspositive columnx=0,y=271fromModel1asaone-dimensionalatmosphere
aslongastheiterationisstartedwithpositivepopulations.This andthehydrogenmodelatom.
is however not true for multigrid iterations. The reason is that
coarsegriditerationshaveamodifiedright-handsideintherate
butthenturnsnegativeuntiliteration22,afterwhichitispositive
equations:
again. For ν = 20 the behaviour is similar, but the number of
2
(cid:40)fk, ifk=(cid:96) iterationswherethepopulationsarenegativeisshorter.
k(nk)= For comparison, in the fourth row we show the relative
W k (fk+1 k+1(nk+1))+ k( k (nk+1)), ifk<(cid:96)
Rk+1 −W W Rk+1 change at the coarsest grid with ν2 = 2 but where we perform
(21) standardMALIiterations(Eq.4)insteadofcoarse-griditerations
(Eq.10).Herethepopulationsarealwayspositive,demonstrat-
Theequationfork < (cid:96)caninprinciplehaveanegativepopula-
ingagainthatMALIiterationskeepsapositivesolutionpositive.
tionasasolution,orproducenegativeestimatesofthesolution
Ifnegativepopulationsoccurafterapplicationofthecorrec-
during iteration. Negative populations on the coarse grids will
tion(Eq.11),wesimplyreplaceitwith10%oftheuncorrected
propagatetothefinestgridanddestroytheconvergence.
populationandcontinueasusualwithpost-smoothingiterations.
By experimentation we found that negative populations at
Ifthishappensatisolatedpointsintheatmospherethenthecon-
the coarse grid typically occur when the error (n n )/n is
notsmoothenough,butwedidnotmanagetogetian−un∞amb∞igu- vergenceofthemultigriditerationisnotaffected,becauseasin-
gleoutlyinggrid-pointrepresentshigh-frequencynoise,whichis
ous definition of what is ’smooth enough’. The error is some-
dampedawayefficiently.However,iftherearetoomanynegative
timesnotsmoothatthelocationofstronggradientsintheatmo-
populationsinthesameareaoftheatmospherethenconvergence
spheric parameters, but sometimes large error fluctuations also
willbeslowerortheiterationmightevendiverge.
occur in smooth parts of the atmosphere models. The error is
typically less smooth in our Hi atom than in the Caii atoms.
This is caused by the larger spectral radius for hydrogen (see 4.11. Initialpopulation
Section 4.13),which lowersthe smoothingcapability of Jacobi
iterations. A good initial approximation for the level populations is criti-
Figure 6 illustrates this in an (cid:96) = 3 test computation using caltoforgoodconvergence.Therearetwopopularmethodsto
a non-smooth 1D atmosphere and the Hi atom. The left-hand initialize the populations: LTE and the zero-radiation-field ap-
columnsshowtheerrorforthefirsttwoenergylevelsoftheatom proximation.Thelattermeanssolvingthestatisticalequilibrium
aftertwoV-cyclesbutwithadifferentnumberofpost-smoothing equations with the radiation field set to zero everywhere in the
iterations.Theerrorforthen=2levelexhibitsastrongspikeat atmosphere.WefoundthatLTE-populationsasaninitialapprox-
a height of 1000 km for ν = 2, but this spike is much smaller imationtypicallyleadstonegativepopulationsatthecoarsegrid.
2
for ν = 10 and ν = 20. In contrast, the error in the n = 1 We suspect that the reason is that the residual and the er-
2 2
populationisrelativelysmooth,andat1000kmheightitisvery rorarelesssmoothusingLTEinitializationthanusingthezero-
small. radiation approximation. To obtain a smooth residual and error
Theblackdotintheleft-handpanelsindicatesthegridpoint manymoreiterationsarerequired.Wedemonstrateourconjec-
that we investigate in more detail in the right-hand panels. The turebycomparingthesmoothingnumber(Hackbush1985;Fabi-
right-handpanelsshowtherelativechangeinthen=1andn=2 aniBendichoetal.1997)asfunctionofMALIiterationsatthe
populationsasafunctionofiterationatthecoarsestgrid(k=1). fine grid for both initializations. The smoothing number is de-
The corrections are small for n = 1, and the n = 1 population finedas
remainspositiveforallcoarse-griditerations,asexpectedfrom
thesmallerrorseenintheleft-handpanels.
ρ (i)=max [n n ] , (22)
L 2 i
Then=2populationbehavesremarkablydifferent.Forν2 = |D − ∞ |
2,thefirstcoarse-griditerationpredictsarelativechangeof-1.7 wheretheoperator isthesecond-orderderivateandiisnum-
2
D
andthusproducesanegativepopulation.Successivecorrections ber of MALI iterations. In Figure 7, we show that the zero-
aresmallerbutthepopulationremainsnegative.Whenν = 10 radiationinitializationleadstoamuchsmoothererrorthanLTE
2
(secondrow),thepopulationispositiveforthefirst3iterations, initialization. A similar result is obtained for the residual. We
Articlenumber,page7of12
A&Aproofs:manuscriptno.accepted_aa
Hin=1 0.2 Hin=2 ×10−7 Hin=1 Hin=2
1.2 4 Postive 10 Postive
∞1.0 0.0 Negative 8 Negative
/n 0.8 0.2 2 6
ν2 =2 )n−∞00..46 −−0.4 /δnn 20 24
(ni 0.2 −0.6 −4 0
0.0 − 2
0.8 −
−
0.1 ×10−7
1.0 2
0.0 4
∞0.8 0.1 0
ν2 =10 /)nn−∞00..46 −−−00..32 /δnn 02 −42
(ni 00..02 −−00..54 −−42 −−6
0.6 8
− −
0.1 ×10−7
0.8
0.0 4 3
n∞0.6 0.1 2 2
/)∞0.4 −0.2 /n 0 1
ν2 =20 n−0.2 −0.3 δn 2 0
ni − −
( 0.0 −0.4 −4 −1
0.5 2
− −
1000 2000 3000 1000 2000 3000 2 ×10−7
Height[Km] Height[Km]
1
10
0
n 1
RHSasEq.5 /n −2 5
δ −
3
−
−4 0
5
−
0 5 10 15 20 25 30 0 5 10 15 20 25 30
ν3 ν3
Fig. 6. Dependence on the number of post-smoothing iterations of the occurrence of negative populations at the coarsest grid. The left-hand
columnsshowtherelativeerroratthefinestgridaftertwofullV-cycles.Theright-handcolumnsshowtherelativechangeofthepopulationatthe
coarsestgridinasinglegridpointintheatmosphere.Reddotsindicatepositivepopulations;bluedotsnegativepopulations.Thechosengridpoint
isindicatedbytheblackdotintheleft-handcolumns.Foreachrowofpanelsadifferentnumberofpost-smoothingiterationsisapplied.Firstrow:
ν =2;secondrow:ν =10;thirdrow:ν =20.Inthefourthrowweshowtherelativechangeofthepopulationatthecoarsestgridwithν =2,
2 2 2 2
butwesolveEq.4insteadofEq.10,i.e.,weperformstandardMALIiterationsinsteadofcoarse-griditerations.Thecomputationwasperformed
ina1Dplane-parallelatmosphereconstructedfromcolumn x = 0,y = 244ofModel1,withtheHiatom.Theothermultigridparametersare
ν =2,ν =32,full-weightingrestrictionandlinearinterpolation.
1 3
therefore use zero-radiation initialization for all our computa-
tions.
S =t /ti , (23)
i MALI MG
4.12. Selectionofpreandpostsmoothingparameters whereti isthecomputingtimerequiredtoperformiV-cycles
MG
(withanassociatederrorE ),andt thecomputingtimeto
i, MALI
reach the same error. A sma∞ll spectral radius means one needs
Theconvergenceandstabilityofthemultigridmethoddepends
fewerV-cyclestoreachconvergence,butthisdoesnotnecessar-
on the selection of the multigrid parameters ν ,ν and ν (see
1 2 3 ily mean that S also increases because one might need many
Sec.2.2).ToensureagoodinitialsmoothingbeforethefirstV- i
iterations (as defined by the value of ν ,ν and ν ) within a V-
cycle,wedefineanadditionalparameter,ν ,whichisthenum- 1 2 3
0 cycletoobtainthesmallspectralradius.
berofMALIiterationsperformedafterinitialization.Thismeans
We found that every atmosphere and model atom behaves
that at the finest grid level we perform ν iterations during the
0 differently,anditisimpossibletogiveasetofparametersthatal-
firstV-cycle,butuseν duringalllaterV-cycles.
1 waysprovidesthefasteststablemultigriditeration.Empirically
Inordertoanalyzethemultigridbehaviorasfunctionofthe we found that two to four pre-smoothings (ν = 2–4) worked
1
νparametersweusetwometrics.Thefirstisthespectralradius. wellinallourtestcases,butforModel1and2withthehydro-
Theothermetricisthespeed-upS ,whichwedefineastheratio gen atom we required extra initial smoothing for the first cycle
i
Articlenumber,page8of12
JohanP.BjørgenandJorritLeenaarts:Numericalnon-LTE3Dradiativetransferusingamultigridmethod
1.0
2.0 ν2=2
ν2=3
1.8 )ρsr 0.8 νν22==45
/tLIMG 1.6 radius( 0.6 ννν222===678
tMA 1.4 ctral 0.4
e
1.2 Sp 0.2
1.0 0.0
1 2 3 4 5 6 7 1 2 3 4 5 6 7
V-cycle V-cycle
Fig.8.Speed-up(left-handpanel)andspectralradius(right-handpanel)
ofmultigridasfunctionofthenumberofV-cycles.Thecomputations
wereperformedusingcolumn x = 0,y = 271fromModel1asaone-
dimensionalatmosphereandthethree-levelCaiimodelatom.Thedif-
Fig.12. TheconvergencepropertiesforatmosphereModel2withsix-
ferent colors indicate a different number of post-smoothing iterations
levelCaiiasfunctionofV-cycle.TheE iscalculatedbasedonthe
ν . The other multigrid parameters are: three grids, ν = 2, ν = 32, approx
2 1 3 relationbetweenthespectralradiusandmaximumrelativepopulation
linearinterpolationandfull-weightingrestriction.
change(Eq.18).Thecomputationwasperformedwiththreegrids,full-
weighting,linearinterpolationusingν =20,ν =4,ν =8,ν =32.
0 1 2 3
(ν 15).Atthecoarsestgrid,wefoundthatν 32isagood
0 3
≥ ≈
number, independent of the problem. Iterating further did not
radiusof0.1-0.2(Trottenbergetal.2000).Weobtainedaspectral
significantlyinfluencetheconvergencespeedorstability.
radius in the range of 0.3-0.4 with the pre and post-smoothing
Themostsensitiveparameteristhepost-smoothingnumber.
parametersasdiscussedinSec4.12usingEq.17.
OurtestswiththeCaiiatomsinModel1andFAL-Cshowedsta-
The spectral radius of one-grid MALI iterations increases
bleandfastconvergencewithν = 2,butforModel2weused
2 withdecreasinggridspacing(Olsonetal.1986).Combinedwith
ν =8toavoidexcessiveoccurrenceofnegativepopulations.In
2 Eq.18,thismeansthatoneneedstoiteratelongerandlongerto
Fig.8weshowthespeedupandspectralradiusina1Dcalcula-
reach the same level of accuracy as the grid spacing gets finer.
tion where we used a column from Model 1 as a plane-parallel
Incontrast,thespectralradiusofmultigridremainsconstantir-
1Datmosphere.Thespectralradiusdecreaseswithincreasingν ,
2 respective of the grid spacing (Steiner 1991; Fabiani Bendicho
butthespeedupgetssmallerowingtotheincreaseincomputa-
et al. 1997; Trottenberg et al. 2000), so that the speed-up of
tionalcostperV-cycle.
multigridisexpectedtoincreasewithdecreasinggridspacing.
For hydrogen we found that one needs many more post
The relation between the maximum relative correction and
smoothing iterations (ν = 25). If a smaller number is chosen
2 the error (Eq. 18) is valid for one-grid MALI iteration (Auer
then the iteration at the coarsest grid tends to produce large ar-
etal.1994).FabianiBendichoetal.(1997)showedthatitholds
easwithnegativepopulations.Wefoundthatthisisrelatedtothe
formultigridinasmooth2Datmospheretoo.The3Dmultigrid
smoothnessoftheerroratthefinestgrid.InFigure 6weillus-
computations presented in Šteˇpán & Trujillo Bueno (2013) did
tratethis.Thetwoleft-handcolumnsshowtheerrorforthefirst
not test the validity of the relation in 3D. Therefore we tested
twolevelsofourhydrogenatomatthestartofthethirdV-cycle
the validity in atmosphere Model 1. The results are shown in
for computations with a different number of post-smoothings
Figure 9, where we compare the true error, the error estimate
(ν =2,10and20).Theerrorgetssmootherwithincreasingν .
2 2 fromEq.18andthemaximumrelativecorrection.For Caiiwe
The right-hand columns show the sign of the level populations
findthattheformulastillholds.Wealsofindthatthemaximum
andthemaximumcorrectionasfunctionoftheiterationnumber
relative correction is equal to the error. This is accidental, be-
atthecoarsestgridduringthethirdV-cycle.Forν = 2anega-
2 cause the spectral radius is close to 0.5. For different pre-and
tivepopulationoccursalreadyduringthefirstcoarse-griditera-
post-smoothingparametersadifferentresultwouldbeobtained.
tionandthisremainssoforallotheriterations.Forν =10neg-
2 Also note that the maximum relative correction in Eq. 18 must
ative populations occur between iteration 3 and 17 after which
beforthewholeV-cycle,andnotfortheindividualMALIiter-
all populations become positive. For ν = 20 negative popula-
2 ations in the post-smoothing at the finest grid. In the hydrogen
tions occur only during iteration 4 to 12. For comparison we
casetheapproximationunderestimatesthetrueerrorbyafactor
showinthefourthrowiterationsatthecoarsestgridforν = 2
2 2. Again the maximum relative correction is just similar to the
wherewesolveEq.4insteadofEq.10,i.e.,weperformstandard
errorbecauseofthespectralradiusbeingcloseto0.5.TheCaii
MALIiterationsinsteadofcoarse-griditerations.Herenonega-
computation reaches the linear regime after five V-cycles. For
tive populations occur, demonstrating that negative populations
Hiitrequiresonlyonecycle.Thisisowingtothemuchhigher
occur owing to the right-hand-side of Eq. 10 and not owing to
number of post-smoothing iterations for hydrogen: ν = 25
thecoarsegridspacing. 2,HI
andν =2.
2,CaII
4.13. Spectralradiusandconvergencecriteria
5. Results
Thespectralradiusforthethree-levelCaiiwithMALIone-grid
iterationinModel1isρ = 0.971,forandHiatomsitisρ = In Table 1 we show a summary of our 3D non-LTE multigrid
sr sr
0.9923.Anoptimalmultigridmethodshouldachieveaspectral computations.Weobtainedaspeed-upof3.3forCaiiand4.5for
Articlenumber,page9of12
A&Aproofs:manuscriptno.accepted_aa
Fig.9. ConvergencebehaviorofmultigridshowingthemaximumrelativeerrorandthemaximumrelativeerrorestimatedbyEq.18foreach
V-cycle,aswellasthemaximumrelativechangeinpopulationperV-cycle.Theupperpanelsshowsthethree-levelCaiiatomandlowerpanels
showsHi,bothforAtmospheremodel1.Forthree-levelCaiiweusedν = 2,ν = 2,ν = 32andforHiweusedν = 15,ν = 2,ν = 25,
1 2 3 0 1 2
ν =32.
3
Table1.Summaryof3Dnon-LTEmultigridruns.
Atmosphere Atom Method ν ,ν ,ν ,ν ,ν Speed-up
FMG 0 1 2 3
Model1 3-levelCaii MALI 1
3-gridMG –,0,2,2,32 3.3
3-gridfull-MG 16,0,2,2,32 6
6-levelHi MALI 1
3-gridMG –,15,2,25,32 4.5
Model2 6-levelCaii 3-gridMG –,20,4,8,32 unknown
HiwithatmosphereModel1andstandardmultigrid.Usingfull tainingareferencesolutionforthismodelwastoocomputation-
multigrid we obtain a speed-up of 6 for Caii, an improvement ally expensive. Therefore, we estimate the convergence based
ofafactor1.8comparedtoordinarymultigrid.Asimilarresult on Eq. 18, which we showed to be accurate for Caii in Model
for a smooth atmosphere was found by Fabiani Bendicho et al. 1 (see Fig. 9). The computation with Model 2 proves that 3D
(1997). In the hydrogen case we needed 25 post-smoothing it- non-LTE radiative transfer using multigrid can handle different
erationstoavoidtoomanyoccurrencesofnegativepopulations. atmospheresandgridspacings.
This considerably reduced the speed-up. The convergence be-
haviorofthesethreerunsisshowninFigure11.
Figure 10 shows the brightness temperature in the Caii K
We also performed multigrid computation in atmosphere linecoreforModel1and2.TheimagecomputedfromModel2
Model2,totestwhethermultigridcanhandleafinergridspac- showssubstantiallymorefine-structure.Thishighlightstheneed
ing(∆ =32km,∆ =13km),andadifferentatmosphere.Ob- forhigher-resolutionradiation-MHDsimulationsandthecapac-
x,y z
Articlenumber,page10of12