BLOCK-STRUCTURED ADAPTIVE MESH REFINEMENT ON - AIP

20 Pages·2012·2.76 MB·English

SIAMJ.SCI.COMPUT. (cid:2)c 2012SocietyforIndustrialandAppliedMathematics
Vol.34,No.3,pp.C102–C121
BLOCK-STRUCTURED ADAPTIVE MESH REFINEMENT ON
CURVILINEAR-ORTHOGONAL GRIDS∗
UDO ZIEGLER†
Abstract. An adaptive mesh refinement scheme for the equations of magnetohydrodynamics
oncurvilinear-orthogonalgridsispresented. Ablock-structuredrefinementapproachisusedsharing
important properties with the hybrid central-upwind/constrained-transport single grid solver like
angularmomentumconservationoncylindricalandsphericalgridsandthepreservationofmagnetic
field solenoidality. The adaptive mesh framework is part of the astrophysical fluid dynamics code
NIRVANAwhichiscontinuouslydevelopedandmaintainedbytheauthor. Thecodeisparallelizedon
the basis of the MPI library, and dynamic load balancing inadaptive meshsimulations is achieved
by space-filling curve domain decomposition techniques. For various two- and three-dimensional
magnetohydrodynamic problemsadaptivemeshfunctionalityandefficiencyaredemonstrated.
Key words. adaptivemeshrefinement,magnetohydrodynamics, numerics: centralschemes
AMS subject classifications. 65M08,65M50
DOI. 10.1137/110843940
1. Introduction. Adaptive mesh refinement (AMR) is a choice method for an
efficient spatial discretization of multiscale computational problems. Among the va-
rieties of existing refinement strategies a block-structured (or patch-structured) ap-
proachhasshowntobeparticularlypopularforthenumericalsolutionoftheequations
ofhydrodynamicsandmagnetohydrodynamics(MHD). Inblock-structuredAMRnu-
merical cells are not individually divided into finer cells, but finer cells are created in
terms of logically rectangular patches which are overlaid the coarser cells. This way
a hierarchy of refinement levels made up of patches can be established in order to
meetresolutionrequirements[3], [5], [4]. Inrecentyearsa couple ofMHD codeshave
been released which follow the block-structured refinement methodology. To name a
few, these arethe FLASH code [10], the SFUMATO code[22], the PLUTOcode [23],
the ENZO code [25], and the AMRVAC code [15], [13]. Some codes like ENZO and
AMRVAC have their own AMR implementation. Others are built on external AMR
packages. For instance, the FLASH code uses the PARAMESH tool [21], whereas
AMR in PLUTO relies on the Chombo library [8]. The use of self-similar blocks
having a fixed number of grid points and a fixed refinement ratio of 2 is a famous
realization of block-structured AMR resulting in a block-tree data structure. Such a
strategyis exploitedby AMRVAC as well as in the astrophysicalfluid dynamics code
NIRVANA presented in this work.
Mostblock-structuredAMRapplicationspublishedintheliteratureareinCarte-
siangeometryusingCartesiangrids. Applicationsbasedontheuseofcylindricalgrids
or spherical grids are, by far, less frequently found especially for three-dimensional
(3D) nonaxisymmetric problems. One reason might be because coding of algorithms
suffersfromanincreasedcomplexityassociatedwiththosegeometries. Afurthercom-
plication represents the intrinsic grid anisotropy caused by the converging grid lines
∗Submitted to thejournal’sSoftware andHigh-PerformanceComputing section August9, 2011;
acceptedforpublication(inrevisedform)February23,2012;publishedelectronicallyJune12,2012.
http://www.siam.org/journals/sisc/34-3/84394.html
†Leibniz-Institutfu¨rAstrophysikPotsdam,D-14482Potsdam,Germany([email protected]).
C102
ADAPTIVEMESHREFINEMENTONCURVILINEARGRIDS C103
at the coordinate singularities degrading solution efficiency and accuracy already on
an uniform grid.
Important properties of the MHD equations which are fulfilled in discrete form
by the numerical scheme on a single grid must be guaranteed on an adaptive grid
as well. Critical issues in this respect requiring extra care are the prolongation and
restrictiontechniquesforsolutiontransferfromcoarse-to-finegridsandfine-to-coarse
grids, and concern mesh synchronization tasks at coarse-fine grid interfaces. In the
numerics of MHD equations the solenoidality constraint for the magnetic field, in
addition, needs special attention. Parallelism of the AMR scheme is a further issue
involving enhanced algorithmic complexity.
Inthispaperablock-structuredmeshrefinementstrategyforcurvilinear-orthogo-
nalgridsisdiscussedasithasbeenrecentlyrealizedintheNIRVANAcode[24]. NIR-
VANA is a numerical code for the approximate solution of the time-dependent, com-
pressible MHD equations. In orthogonal-curvilinear coordinates (x,y,z) the metric
tensor has diagonal form, H = diag(h2,h2,h2). For cylindrical coordinates (z,R,φ)
x y z
and spherical coordinates (r,θ,φ) the metric scale factors are given by (hx,hy,hz)=
(1,1,R) and by (hx,hy,hz)= (1,r,rsinθ), respectively. In the trivial case of Carte-
siancoordinates the scale factors are (hx,hy,hz)=(1,1,1). The governingequations
are
1
∂tu=− [∂x(hyhzfx)+∂y(hzfy)+∂z(hyfz)]+S,
hyhz
1
∂tBx =− [∂y(hzEz)−∂z(hyEy)] ,
hyhz
1
∂tBy =− [∂zEx−∂x(hzEz)] ,
hz
1
(1.1) ∂tBz =− [∂x(hyEy)−∂yEx]
hy
with flux functions fx,fy,fz,
⎛ ⎞
(cid:5)vx
⎜ ⎟
⎜⎜ (cid:5)vx2+pt−Bx2/μ ⎟⎟
fx(u,B)=⎜⎜ (cid:5)vxvy −BxBy/μ ⎟⎟ ,
⎜ ⎟
⎝ hz((cid:5)vxvz −BxBz/μ) ⎠
(e+pt)vx−Bx(v·B)/μ
⎛ ⎞
(cid:5)vy
⎜ ⎟
⎜⎜ (cid:5)vxvy −BxBy/μ ⎟⎟
fy(u,B)=⎜⎜ (cid:5)vy2+pt−By2/μ ⎟⎟ ,
⎜ ⎟
⎝ hz((cid:5)vyvz −ByBz/μ) ⎠
(e+pt)vy −By(v·B)/μ
⎛ ⎞
(cid:5)vz
⎜ ⎟
⎜⎜ (cid:5)vxvz −BxBz/μ ⎟⎟
fz(u,B)=⎜⎜ (cid:5)vyvz −ByBz/μ ⎟⎟
⎜ ⎟
⎝ hz((cid:5)vz2+pt−Bz2/μ) ⎠
(e+pt)vz −Bz(v·B)/μ
C104 UDOZIEGLER
and geometric source term S:
⎛ ⎞
0
⎜ (cid:8) (cid:9) ⎟
⎜ 1 ∂hy (cid:5)(v2+v2)+2p+B2/μ ⎟
⎜⎜ hy (cid:10)∂x y z x (cid:11) ⎟⎟
S(u,B)=⎜⎜⎜ h1y h1z ∂∂hyz((cid:5)vz2+pt−Bz2/μ)− ∂∂hxy((cid:5)vxvy−BxBy/μ) ⎟⎟⎟ .
⎝ 0 ⎠
0
Here, u = ((cid:5),(cid:5)vx,(cid:5)vy,hz(cid:5)vz,e) is the vector of primary hydrodynamical variables
with (cid:5) the mass density, v = (vx,vy,vz) the fluid velocity, and e the total energy
density. The magnetic field is denoted by B = (Bx,By,Bz), and the electric field E
is defined by E = −v×B. The total pressure pt = p+B2/2μ is the sum of gas
pressure p and magnetic pressure with μ the mag(cid:12)netic permeability. Th(cid:13)e equations
are complemented by an ideal gas law p=(γ−1) e−(cid:5)|v|2/2−|B|2/2μ with adia-
batic coefficient γ and by the solenoidality property of the magnetic field, ∇·B = 0.
Note that the fourth component of the u-equation (for u = hz(cid:5)vz) has no source
term expressing the conservation of angular momentum in cylindrical and spherical
coordinates.
The organization of the paper is as follows. First, in section 2 NIRVANA’s
MHD solver is summed up. This is followed by the description of NIRVANA’s AMR
paradigm in section 3. A verification of the AMR scheme is then given in section 4.
Severalnumericalexamplesintwo-dimensional(2D)and3Dcylindricalandspherical
geometry are presented. Finally, conclusions are drawn in section 5.
2. Single grid solver: Summary. In block-structuredAMR updating the nu-
merical solution from a time level tn to a next time level tn+1 =tn+δt (where δt is
the discrete time-step) consists of applying a single grid solver to each grid patch in
the mesh hierarchy separately and then matching those individual solutions in space
and, eventually, in time in case different time-steps are used for patch updates.
In the following NIRVANA’s single grid MHD solver is reviewed insofar as nec-
essary for an understanding of the applied AMR philosophy. The very details of the
solverarenotrepeatedhere, andthe readerisreferredtoZiegler[33] foramorecom-
prehensive discussion. In NIRVANA the set of equations (1.1) is solved by a hybrid
method consisting of a finite-volume (FV) approximation of Godunov type for the
u-update (Euler subsystem including the Lorentz force) and a constrained-transport
(CT)schemefortheB-update(inductionequation). Byconstruction,theCTscheme
guaranteesthatthe time derivativeofmagneticdivergencevanishes preciselyina FV
sense [30], [20], [2], i.e., d∇·B/dt = 0, where the overbar denotes volume averaging
overagridcell. TheintegralformofMaxwell’ssolenoidalityconditionisthereforeau-
tomaticallyfulfilledcellwiseifsocellwiseinitially. Additionalgriddivergencecleaning
procedures (see, e.g., [9], [6]) are redundant in CT schemes.
The FV form of (1.1) reads componentwise as
u˙i,j,k(t)=−[δAx.,j,kFx.,j,k]iipm +[δAyi,.,kFyi,.,k]jjpm +[δAzi,j,.Fzi,j,.]kkpm +Si,j,k
δVi,j,k
(2.1)
and describes the time rate of change of volume-averages
(cid:14)
1
ui,j,k = u(x,t)dV
δVi,j,k C
i,j,k
ADAPTIVEMESHREFINEMENTONCURVILINEARGRIDS C105
in cells Ci,j,k = [xim,xip] × [yjm,yjp] × [zkm,zkp] of a uniform grid represented by
running indices i,j, and k and grid spacings δx = xip −xim, δy = yjp −yjm, and
δz =zkp−zkm. Halfindicesim =i−1/2,ip =i+1/2andanalogindicesfortheother
directions mark cell nodes with coordinates(xim,yjm,zkm). The geometriccell center
has coordinates (xi,yj,zk). The cell volume is denoted by δVi,j,k. The areas of the
cell bounding x-faces aredenotedby δAxim,j,k forthe cell faceatxim andδAxip,j,k for
the cellface atxip, respectively. Analogdefinitions applyforothercellfaces. Fx, Fy,
andFz arethe face-averagedflux functions. Asubscriptedsquarebracketdenotesan
index difference; e.g., [δAx.,j,kFx.,j,k]iipm = δAxip,j,kFxip,j,k −δAxim,j,kFxim,j,k means
the flux difference betweenthe upper andlowerx-faces. The numerical fluxes areap-
proximatedfromtheface-averagedfluxfunctionsby(i)computingthefacialintegrals
with the help of a cubature formula of degree 2 and by (ii) utilizing central-upwind
techniques for stability. Essentially, the method is an extension of the second-order
version of the semidiscrete scheme developed in [18] to curvilinear-orthogonal grids
(see also the work in [14], [17]). The resulting flux formulae are equivalent to HLL
fluxes (see equation (12) in [33]).
Point values of variables u needed in the flux computation process are obtained
from reconstructing the solution based on the known set of averages {ui,j,k}. The
reconstruction function u˜(x) in a cell Ci,j,k is given by
(2.2) u˜(x)=ui,j,k+((cid:2)∂xu)i,j,k(x−xi)+((cid:2)∂yu)i,j,k(y−yj)+(cid:2)(∂zu)i,j,k(z−zk),
where xi, yj are cell-volume-averaged coordinates. Explicit expressions for those in
cylindrical and spherical coordinates are given in Table 1 in [33]. The numerical
(cid:2) (cid:2) (cid:2)
derivatives (∂xu)i,j,k, (∂yu)i,j,k, (∂zu)i,j,k are approximations to the exact partial
derivativessubjecttoalimiterfunctionforoscillationcontrolnearflowdiscontinuities.
The harmonic mean limiter of [29] has been used in all of the test problems in this
paper. The reconstruction is conservative,
(cid:14)
ui,j,kδVi,j,k = u˜(x)dV ,
C
i,j,k
and has the obvious property ui,j,k =u˜(xi,yj,zk); i.e., the cell-averageis identical to
the point value at coordinates (xi, yj,zk).
The CT scheme for the (staggered) magnetic field components reads as
(cid:15) (cid:16)
(cid:3)B˙x(cid:4)xim,j,k(t)= δAxi1m,j,k (cid:15)−δz[hzim,.Ezim,.,k]jjpm +δyhyim[Eyim,j,(cid:16).]kkpm ,
(cid:3)B˙y(cid:4)yi,jm,k(t)= δAyi1,jm,k (cid:15)−δx[Exi,jm,.]kkpm +δz[hz.,jmEz.,jm,k](cid:16)iipm ,
(2.3) (cid:3)B˙z(cid:4)zi,j,km(t)= δAzi1,j,km −δy[hy.Ey.,j,km]iipm +δx[Exi,.,km]jjpm
and describes the time rate of change of face-averages
(cid:14)
1
(cid:3)Bx(cid:4)xim,j,k = δAxim,j,k (cid:14)δAxim,j,kBx(x,t)dA,
(cid:3)By(cid:4)yi,jm,k = δAyi1,jm,k (cid:14)δAyi,jm,kBy(x,t)dA,
1
(cid:3)Bz(cid:4)zi,j,km = δAzi,j,km δAzi,j,km Bz(x,t)dA,
C106 UDOZIEGLER
where (cid:3).(cid:4)x(y,z) means face-averaging over δAx(y,z). The divergence-free condition in
discrete form becomes
ddt(cid:17)(cid:8)(cid:3)Bx(cid:4)x.,j,kδAx.,j,k(cid:9)iipm +(cid:10)(cid:3)By(cid:4)yi,.,kδAyi,.,k(cid:11)jjpm +(cid:8)(cid:3)Bz(cid:4)zi,j,.δAzi,j,.(cid:9)kkpm(cid:18)=0.
The edge-averagedelectric field E appearing in equations (2.3) is approximated with
the help of a semidiscrete evolution-projection method originally developed for dis-
cretizingequationsofHamilton–Jacobitype[7], [18]. Theresultingformulaeareeach
equivalent to a composition of one-dimensional (1D) HLL Riemann solvers in the re-
spective 2D subspace reflecting the 2D nature of the individual components of the
induction equation (see equations (20)–(22) in [33]). Also, this approach naturally
reduces to a correct 1D upwinding for grid-aligned flows (see also [11], [12] for an
alternative ansatz).
The reconstruction functions for the magnetic field components are given by
B˜x(x)= 21(cid:15)(cid:3)Bx(cid:4)xip,j,k+(cid:3)Bx(cid:4)xim,j,k(cid:16)+(cid:8)(cid:3)Bx(cid:4)x.,j,k(cid:9)iipm xΔ−xxi
(cid:19)(cid:15) (cid:16) (cid:15) (cid:16) (cid:20)
+ ∂(cid:3)yBx x−xim + ∂(cid:3)yBx xip −x (y−(cid:3)y(cid:4)xj)
ip,j,k Δx im,j,k Δx
(cid:19)(cid:15) (cid:16) (cid:15) (cid:16) (cid:20)
+ ∂(cid:3)zBx x−xim + ∂(cid:3)zBx xip −x (z−zk),
ip,j,k Δx im,j,k Δx
(cid:19)(cid:15) (cid:16) (cid:15) (cid:16) (cid:20)
B˜y(x)= ∂(cid:3)xBy y−yjm + ∂(cid:3)xBy yjp −y (x−(cid:3)x(cid:4)yi)
i,jp,k Δy i,jm,k Δy
(cid:15) (cid:16) (cid:10) (cid:11)
+ 21 (cid:3)By(cid:4)yi,jp,k+(cid:3)By(cid:4)yi,jm,k + (cid:3)By(cid:4)yi,.,k jjpm yΔ−yyj
(cid:19)(cid:15) (cid:16) (cid:15) (cid:16) (cid:20)
+ ∂(cid:3)zBy y−yjm + ∂(cid:3)zBy yjp −y (z−zk),
i,jp,k Δy i,jm,k Δy
(cid:19)(cid:15) (cid:16) (cid:15) (cid:16) (cid:20)
B˜z(x)= ∂(cid:3)xBz z−zkm + ∂(cid:3)xBz zkp −z (x−(cid:3)x(cid:4)zi)
i,j,kp Δz i,j,km Δz
(cid:19)(cid:15) (cid:16) (cid:15) (cid:16) (cid:20)
+ ∂(cid:3)yBz z−zkm + ∂(cid:3)yBz zkp −z (y−(cid:3)y(cid:4)zj)
i,j,kp Δz i,j,km Δz
(2.4) + 21(cid:15)(cid:3)Bz(cid:4)zi,j,kp +(cid:3)Bz(cid:4)zi,j,km(cid:16)+(cid:8)(cid:3)Bz(cid:4)zi,j,.(cid:9)kimp zΔ−zzk ,
(cid:3)
where ∂yBx and other cross derivative terms are slope-limited approximations to
the exact derivatives. (cid:3)x(cid:4)y ((cid:3)x(cid:4)z,(cid:3)y(cid:4)x,(cid:3)y(cid:4)z) is the face-averaged x(x,y,y)-coordinate
i i j j
where averaging is over the cell face δAy (δAz,δAx,δAz). As can be shown by explicit
ADAPTIVEMESHREFINEMENTONCURVILINEARGRIDS C107
calculation this B-reconstructionis flux-conserving in the sense that
(cid:14)
1
δAxim,j,k (cid:14)δAxim,j,kB˜x(x)dA=(cid:3)Bx(cid:4)xim,j,k,
δAyi1,jm,k (cid:14)δAyi,jm,kB˜y(x)dA=(cid:3)By(cid:4)yi,jm,k,
1
δAzi,j,km δAzi,j,km B˜z(x)dA=(cid:3)Bz(cid:4)zi,j,km
and is consistent with the solenoidality property of the CT scheme,
δVi,j,k(∇·B˜)i,j,k =(cid:8)(cid:3)Bx(cid:4)x.,j,kδAx.,j,k(cid:9)iipm +(cid:10)(cid:3)By(cid:4)yi,.,kδAyi,.,k(cid:11)jjpm +(cid:8)(cid:3)Bz(cid:4)zi,j,.δAzi,j,.(cid:9)kkpm .
The semidiscrete scheme (2.1) and (2.3) is finally solved using a time-explicit
Runge–Kutta integrator. NIRVANA has the choice between the standard second-
orderaccurate Runge–Kutta method andthe more stable third-ordermethod of [27].
3. AMRoncurvilinear-orthogonalgrids. Amajortechnicaltaskanyblock-
structured AMR implementation has to master concerns the managementof the grid
hierarchy which consists of a set of individual grid blocks with different resolutions.
For efficiency reasons it is inevitable that such a block management system also pro-
videsinformationontheconnectivityofablockwithneighboringblocks,childblocks,
andparentblocks. Further problemsonehasto dealwith areassociatedwiththe FV
framework for the MHD system of conservation laws. Important issues are (i) how
to evolve the grid hierarchy one time-step within the Runge–Kutta multistage inte-
grationprocess, (ii) how to match solutions at coarse-finegridinterfaces in a manner
consistentwiththeFVapproach,and(iii)howtotransferthecoarsegridsolutiontoa
generatedfiner grid. Intwo previouspapersby Ziegler[31], [32] all suchaspectshave
already been discussed for AMR on Cartesian grids. Because they build upon the
same principles, majorparts ofthe CartesianAMR implementationcarryoverto the
more general case of AMR on curvilinear-orthogonal grids. A summary of the basic
ideas is givenbelow. However,changesto the algorithmbecome necessaryinthe flux
and electric field synchronization process at coarse-fine grid interfaces as well as in
therestrictionandprolongationoperationsofsolutiontransfer. Infact,theprolonga-
tion algorithm for the magnetic field used in Cartesian geometry failed in preserving
the solenoidality property on general curvilinear-orthogonalgrids and, therefore, has
been replaced by a new hierarchical prolongation approach. The failure is because
the Cartesian algorithm relies on simpler interpolation violating flux conservationon
generalcurvilinear-orthogonalgrids. Suchmodificationsaredescribedinmoredetail.
3.1. General strategy. Block-structured AMR in NIRVANA is based on the
use of generic blocks having a fixed size of four cells per coordinate direction. In the
beginningasetofbaselevelblocksspansthecomputationaldomain. Refinementsare
then achieved by a hierarchical nesting of blocks assuming a fixed refinement ratio
of 2 between adjacent refinement levels in all coordinate directions. A 2D (3D) finer
block alwayscoversa quad(oct) of2×2 (2×2×2)underlying coarserblock cells. A
system of nested blocks on each base level block is established building a block-tree
structure. Block connectivity within such a tree as well as connectivity to spatially
neighboring blocks is efficiently handled by pointers mediating instant data access.
Load balancing in parallel simulations is achieved by an equal distribution of blocks
among threads following space-filling curve techniques [26]. A space-filling curve of
C108 UDOZIEGLER
Hilbert type serves for ordering base level blocks and an adapted curve of Morton
type for the order of blocks within a base block’s tree.
A finer block is generatedaccording to some mesh refinement criteria. The stan-
dard criterion in the NIRVANA code is a convex combination of criteria based on
the relative strengths of first (respectively, second) derivatives of certain variables u.
The expression for Cartesian AMR is detailed in [32]. It needs only slight adapta-
tionforcurvilinearAMR.Accordingtothisstandardcriterion,meshrefinementtakes
place if a prescribed thresholdis exceeded. Those u-thresholdsare usually somewhat
problem-dependent parameters. In addition to the derivatives-based criterion, mesh
refinement in NIRVANA may also be controlled by a Jeans-length-based criterion in
case of self-gravitational flows or a Field-length-based criterion for flows containing
thermally unstable gas. Here, in the test problems the standard criterion is used
with the density, velocity, and magnetic field components being refinementtriggering
variables as the default.
Althoughtheuseofsmall-sized,self-similarblockshasclearadvantagesintermsof
grid adaptivity, mesh organization, and load balancing issues, advancing the solution
onalargenumberofsmallblocksreducescomputationalefficiencyduetotheincreased
synchronizationoverheadatblockinterfaces. Inorderto speedupintegration,blocks
ofa givenrefinementlevelaretherefore temporarilyclusteredto larger-sizedpatches.
In parallel simulations the clustering is done separately for each domain partition.
Ghostcellsareaddedforpatchesinplaceofgenericblocks. Thediscretizedequations
are then solved on these patches. The solution is copied back to generic blocks after
a time-step is complete, and patches are then deallocated. Note that the additional
overhead due to clustering is usually substantially overcompensated by the gain in
computational speed. The adaptive grid is advanced with a single time-step given
by the minimum of CFL-related time-steps computed over the full patch hierarchy;
hence,thereisnorefinementintime. WithinthemultistageRunge–Kuttaintegration
scheme the full patch hierarchy is updated stage-by-stage. Each patch is treated as
a separate initial value problem with given boundary conditions obtained from grid
neighbors (either copied from equal level grids or interpolated from coarser grids) or
defined at physical boundaries. After each Runge–Kutta stage a synchronization of
ghostcellstakesplaceatintralevelpatchinterfacestoretrievesolutionconsistency. In
particular, because of the staggering of magnetic field components those components
atpatchinterfaceslivingondifferentpatches,localorremote,musthaveuniquevalues
without ambiguity. In a restriction operation the less accurate coarse grid solution
is replaced by the fine grid solution in overlapping regions. A flux and electric field
correction step is necessary at fine-coarse grid interfaces to recover conservation and
to maintain the magnetic field solenoidality constraint.
3.2. Prolongation operation. Forcell-averagedquantitiesa conservativepro-
longation of the coarse grid solution onto the cells of a newly generated finer block
canbe constructedentirelyonthe basisofthe reconstructionformula(2.2)forcoarse
cells. For a coarse cell Cc covered by the oct (in three dimensions) of finer cells
i,j,k
Cf with running indices l,m,n, conservation is expressed by
l,m,n
(cid:21)
c c f f
(3.1) u δV = u δV ,
i,j,k i,j,k l,m,n l,m,n
l,m,n
(cid:22)
whereasuperscriptc(f)marksacoarse(fine)cellquantityandδVc = δVf .
i,j,k l,m,n l,m,n
f f
Fine grid cell-averaged values u obeying (3.1) are simply given by u =
l,m,n l,m,n
ADAPTIVEMESHREFINEMENTONCURVILINEARGRIDS C109
u˜c(xl,ym,zn),where(xl,ym,zn)arethecoordinatesofthevolumetriccenterofClf,m,n.
The proof is straightforward:
(cid:21) (cid:21)
f f c f
ul,m,nδVl,m,n = u˜ (xl,ym,zn)δVl,m,n
l,m,n l,m,n
(cid:23) (cid:24)
(cid:21)
c
=uci,j,kδVic,j,k+((cid:2)∂xu)i,j,k xlδVlf,m,n−xiδVic,j,k
l,m,n
(cid:23) (cid:24)
(cid:21)
c
+((cid:2)∂yu)i,j,k ymδVlf,m,n−yjδVic,j,k
l,m,n
(cid:23) (cid:24)
(cid:21)
c
+(cid:2)(∂zu)i,j,k znδVlf,m,n−zkδVic,j,k =uci,j,kδVic,j,k
l,m,n
since the terms in brackets vanish. For example,
(cid:14)
(cid:21) (cid:21)
xlδVlf,m,n−xiδVic,j,k = xlδVlf,m,n− xdV
l,m,n l,m,n Cic,j,k(cid:14)
(cid:21) (cid:21)
= xlδVlf,m,n− xdV
l(cid:21),m,n l(cid:21),m,n Clf,m,n
= xlδVlf,m,n− xlδVlf,m,n =0.
l,m,n l,m,n
The prolongation procedure for the staggered magnetic field components is a bit
more complex. Here, in order to ensure magnetic field solenoidality in each new fine
cell of the oct it is not sufficient to solely take point values from the coarse grid
reconstruction at the proper fine cell positions. However, the formulae (2.4) are used
to compute fine cell values at a coarse cell’s faces wherever there is no neighboring
fine grid from which values are naturally copied. Accordingly, flux-conserving new
fine values at the coarse cell faces are given by the point values of B-reconstruction
at proper positions, i.e.,
(cid:3)Bxf(cid:4)xlm,m,n =B˜xc(xlm,(cid:3)y(cid:4)xm,zn),
(cid:3)Byf(cid:4)yl,mm,n =B˜yc((cid:3)x(cid:4)yl,ymm,zn),
(cid:3)Bzf(cid:4)zl,m,nm =B˜zc((cid:3)x(cid:4)zl,(cid:3)y(cid:4)zm,znm).
The situation is illustrated in Figure 3.1(a). It remains to compute the 12 magnetic
field values at the inner fine cell faces. This is done using a hierarchical prolongation
procedureasfollows. Inafirststepthecoarsecellisvirtuallydividedintotwosubcells
with the cutting plane x = xi through the geometrical cell center (Figure 3.1(b)).
With the assumption that each of the two x-halves is divergence-free in an FV sense
an intermediate (superscript I) face-averaged value (cid:3)BxI(cid:4)xi,j,k at xi can be calculated
from the divergence condition applied to any half. Taking, for instance, the x-lower
half the condition reads as
(cid:3)BI(cid:4)x δAI −(cid:3)Bc(cid:4)x δAc
x i,j,(cid:23)k xi,j,k x im,j,(cid:24)k xim(cid:23),j,k (cid:24)
(cid:21) jp (cid:21) kp
+ (cid:3)Bf(cid:4)y δAf + (cid:3)Bf(cid:4)z δAf =0,
y l,.,n yl,.,n z l,m,. zl,m,.
n jm m km
where the bracketscontain the previously computed fine cell values at the coarsecell
C110 UDOZIEGLER
z y f
Bz
x
f
Bx
(a)
f
By
f
I Bx
(b) Bx
f
By
I
By
(c)
f
Bz
(d)
Fig. 3.1. Magnetic field prolongation strategy: (a) fine cell values at coarse cell faces, (b)
prolongation of interior fine cell x-component, (c) prolongation of interior fine cell y-component,
(d) prolongation of fine cellz-component.
faces. Then, the four fine cell x-components of the magnetic field at xi are given by
(cid:3)Bxf(cid:4)x.,m,n =B˜xI((cid:3)y(cid:4)xm,zn),
where B˜xI(y,z) is a flux-conservative 2D Bx-reconstruction in the plane x= xi. It is
determined from the set {(cid:3)BI(cid:4)x } and, by construction, is equipped with the prop-
x i,j,k (cid:22)
erty (cid:3)BxI(cid:4)xi,j,k =B˜xI((cid:3)y(cid:4)xj,zk). Indeed, it can be shown that m,n(cid:3)Bxf(cid:4)x.,m,nδAfx.,m,n =
ADAPTIVEMESHREFINEMENTONCURVILINEARGRIDS C111
(cid:3)BI(cid:4)x δAI , expressing flux conservation at the intermediate x-face. In the next
x i,j,k xi,j,k
step the coarse cell is further divided with the cutting plane this time at y = yj,
resulting in four subcells (Figure 3.1(c)). Again, applying the divergence condition
to any quarter of the previous halves gives the two intermediate values (cid:3)BI(cid:4)y .
y i±1,j,k
4
Taking, for instance, the y-lower subcell in the x-lower half the condition reads as
(cid:23) (cid:24)
(cid:21) i
(cid:3)Bf(cid:4)x δAf +(cid:3)BI(cid:4)y δAI
x .,m,n x.,m,n y i−1,j,k yi−1,j,k
n (cid:21) im (cid:10) 4 4 (cid:11)
− (cid:3)Bf(cid:4)y δAf + (cid:3)Bf(cid:4)z δAf kp =0,
n y l,mm,n yl,mm,n z l,m,. zl,m,. km
where in (cid:3)BI(cid:4)y the average is taken over the area δAI with half x-range
y i−1,j,k yi−1,j,k
4 4
[xim,xi]. Note that this step needs the previously computed fine values (cid:3)Bxf(cid:4)x.,m,n.
With the help of the computed sets {(cid:3)BI(cid:4)y }, flux-conservative 1D reconstruc-
y i±1,j,k
tion functions B˜yI±(z)along the z-directionan4d throughpoints ((cid:3)x(cid:4)yl,yj) can be con-
structed equipped with the property (cid:3)ByI(cid:4)yi±1,j,k = B˜yI±(zk). Then, the four fine cell
4
y-components of the magnetic field at yj are simply given by taking point values of
the respective reconstruction function:
(cid:3)Byf(cid:4)yl,mm,n =B˜yI±(zn).
Finally, cutting the coarse cell additionally at z = zk the oct of fine cells results
(Figure 3.1(d)). The remaining four fine cell z-components of the magnetic field
at zk are obtained directly by an application of the divergence condition in each of
the four octants in either the z-lower half or z-upper half of the coarse cell. This
completes the prolongation operation which, by construction, is a divergence-free
transfer of the coarse grid magnetic field solution onto the fine grid. Note that the
presented prolongationstrategy is similar to the one described in [19], although both
approachesmaydifferindesigndetails. Acompletelydifferentansatzforadivergence-
free reconstruction is used in the polynomial methods of [28], [1].
3.3. Flux and electric field correction. A conservativerestriction operation
is applied when transferring the more accurate fine grid solution in a hierarchical
manner onto the coarser grid. For cell-averaged quantities the formula is given by
(3.1). As a consequence the computed coarse grid flux and electric field components
lyingdirectlyatcoarse-finegridinterfacesarenolongerconsistentwiththecopiedso-
lution. Inordertoensureconsistencyand,thus,torecoverconservationandmagnetic
field solenoidality a correction step is necessary. Obviously, the coarse grid fluxes at
those interfacesmust be replacedby the sum offine gridfluxes for the quad(in three
dimensions) of fine cell faces covering a coarse cell face:
(cid:21)
Fc −→Fc,corr , δAc Fc,corr = δAf Ff ,
xim,j,k xim,j,k xim,j,k xim,j,k xlm,m,n xlm,m,n
m,n
(cid:21)
Fc −→Fc,corr , δAc Fc,corr = δAf Ff ,
yi,jm,k yi,jm,k yi,jm,k yi,jm,k yl,mm,n yl,mm,n
l,n
(cid:21)
Fc −→Fc,corr , δAc Fc,corr = δAf Ff .
zi,j,km zi,j,km zi,j,km zi,j,km zl,m,nm zl,m,nm
l,m
In analogy, the edge-centered coarse grid electric field components must be replaced
by the sum of fine grid values for the double of fine cell edges covering a coarse cell