David A. Clarke & Jonathan P. Ramsey

AMR
+
ZEUS‐3D
=
AZEuS
 An
Adap/ve
Zone
Eulerian
Scheme


HALIFAX, NOVA SCOTIA, CANADA

David
A.
Clarke
&
Jonathan
P.
Ramsey


Abstract

3. N (n) is the number of coarse (fine) elapsed time steps: tN = tn .

Adaptive Mesh Refinement (AMR) has been incorporated into nearly every major effort in computational astrophysical fluid dynamics, including MHD. The vast majority of these efforts have been with zone-centred and usually Godunov-type schemes, but we show here that AMR can also be adapted to a staggered mesh MHD code such as ZEUS-3D. In this poster, we discuss the challenges a staggered mesh poses to refluxing and interpolation, and how the principles of conservation and the solenoidal condition can be maintained. For our problem in particular, it is critical that all linear and non-linear waves pass through grid boundaries with a minimum of reflection, refraction, and dissipation, and we demonstrate how we have achieved this.

4. The zone volume, face area, edge length, and time step are ∆V , ∆A, ∆x, and ∆t for the coarse grid; δV , δA, δx, and δt for the fine grid.

solves the following system of equations: ∂ρ + ∇ · (ρ #v ) = 0; ∂t ! " ∂#s # # + ∇ · #s #v + (p + pB ) I − B B + Q = 0; ∂t ! " ∂eT # ×B # + #v · Q = 0; + ∇ · (eT + p − pB ) #v − E ∂t # ∂B # = 0, +∇×E ∂t

(1)

5. For illustration, Cartesian coordinates and constant geometric factors (∆V , δA, etc.) are assumed (AZEuS is covariant in xyz, zrφ and rθφ coordinates). Thus, ∆V /δV = $3 , ∆A/δA = $2 , ∆x/δx = $ and ∆t/δt = $ where $ is the refinement ratio, always a power of 2.

(2)

6. The volume integral of a volume-conserved quantity is conserved (e.g.,

(3)

7. The surface integral of an area-conserved quantity is conserved (e.g.,

(4)

8. Discretised flow variables [e.g., P (I, J, K), s1 (i, j, k), etc.] are averages of the actual quantity (ρ, s1 ) over their region of influence (RI). For example:

ρdV = M).

& · dA & = ΦB ). B

- The RI of a 1-face-centred volume-conserved quantity [S1 (I, J, K), Fig. 2b] is the right half of zone (I −1, J, K) plus the left half of zone (I, J, K).

p = (γ − 1)e is the thermal pressure, e is the internal energy density; pB = 12 B 2 is the magnetic pressure (energy density) in units where µ0 = 1;

- The RI of a 3-face-centred area-conserved quantity [B3 (I, J, K), Fig. 2c] is the (I, J, K) 3-face.

I is the identity tensor;

(Berger, M. J. & Colella, P., 1989, J. Comput. Phys., 82, 64.)

A

V

- The RI of a zone-centred volume-conserved quantity [e.g., ρ(i, j, k), Fig. 2a] is zone (i, j, k).

where: #s = ρ#v is the momentum density;

Adaptive Mesh Refinement

!

!

Q is the diagonal artificial viscous tensor where Qii = ρ∆vi (qq ∆vi − ql cs ), qq = 0 if ∆vi > 0;

Adaptive Mesh Refinement (AMR) allows a fluid simulation to be performed with sufficient resolution to follow physically critical phenomena without over-resolving regions (and thus wasting computing resources) that can be followed at lower resolution. It is essentially a bookkeeping algorithm constrained by physics that manages multiple nested and/or overlapping grids of varying resolution in such a way that conserved variables remain conserved. AMR also establishes when a grid must be added or may be removed and effects such modifications.

The restriction step

eT = e + 12 ρv 2 + pB is the total energy density; # = −#v × B # is the induced electric field, E

Restriction depends on how the variable is centred and conserved. For volume-conserved quantities, integrate over the volume of the coarse RI to get:

with all other variables having their usual meaning. ZEUS-3D uses an operator split algorithm for the momentum, a conservative operator un-split algorithm for the total energy, and is covariant for xyz, zrφ, and rθφ coordinates. It is upwinded in entropy and shear Alfv´en waves with (magneto)sonic waves stabilised using a von Neumann artificial viscosity. Interpolation is time-centred and can be performed using piecewise linear or parabolic reconstruction.

AMR accomplishes all this by performing five principle tasks (the “five Res”): 1. Restriction: coarse grid variables are replaced with conservative averages of those in an overlying fine grid;

Our problem...

3. Reconstruction (a.k.a., prolongation): a new fine grid is created from an underlying coarse grid with a suitable locally conservative interpolation scheme;

A proto-stellar jet is launched at the sub-AU scale and propagates to the sub-parsec scale where it is observable (background image). As the jet expands over five orders of magnitude in scale in a roughly self-similar fashion, there is no need to resolve the jet at sub-parsec scales with the same resolution as at sub-AU scales. This suggests fixed nested grids through which the jet propagates and between which all linear and non-linear MHD waves must pass cleanly with minimal reflection, refraction, and degradation. # ZEUS-3D is built on a staggered mesh, in which scalars (ρ, eT ) are zone-centred, primary vectors (#s, B) # J) # are edge-centred (Fig. 1). Conversely, the original AMR is are face-centred, and derived vectors (E, designed for fully zone-centred schemes.

5. Grid restructuring: Based on some physical and/or efficiency criteria, grids are created, eliminated, split, or merged. This step is not considered in this poster.

QdV =

∆V

qdV.

Thus, restriction of a zone-centred quantity (e.g., ρ, eT ; Fig. 2a) is given by: (5)

where $ is the refinement ratio. The restriction of a 1-face-centred quantity (e.g., s1 ; Fig. 2b) is given by: !/2 # !−1 # 1 # !−1 S1 (I, J, K) = 3 ηα s1 (i+α, j +β, k+γ), $ ! β=0 γ=0 α=−

(6)

2

where ηα = 12 for α = ± 2! , 1 otherwise, since only half of the left (i − !2 ) and right (i + 2! ) fine RIs of s1 cover the coarse RI of S1 . “Skin values” of S1 (those penetrating the outer-boundary of the fine grid) are not restricted since only half of their RI is covered by fine zones. Instead, these are refluxed (e.g., Fig. 3). For area-conserved quantities such as B3 , integrate over the area of the coarse RI to get:

The need to pass all waves across grid boundaries and to accommodate a staggered mesh requires significant modifications to the original AMR scheme.

ZEUS-3D

∆V

"

# !−1 # !−1 # 1 !−1 q(i+α, j +β, k+γ), Q(I, J, K) = 3 $ α=0 β=0 γ=0

2. Refluxing: fluxes computed from the coarse grid into a coarse zone adjacent to (but not underlying) a fine grid are replaced by fluxes computed from the fine grid;

4. Boundary restoration: Boundary zones of a fine grid are restored from underlying coarse active zones with a suitable locally conservative interpolation scheme (similar to reconstruction but sufficiently different to be treated as a separate step);

"

"

∆A

B3 dA =

"

∆A

b3 dA.

Thus, restriction of B3 (Fig. 2c) is given by:

(Clarke, D. A., 1996, ApJ, 457, 291; http://www.ica.smu.ca/zeus3d)

Conventions and definitions

# !−1 # 1 !−1 B3 (I, J, K) = 2 b3 (i+α, j +β, k), $ α=0 β=0

ZEUS-3D is the well-known astrophysical fluids solver which, in its ideal MHD mode without self-gravity, 1. Coarse- and fine-zone quantities are indicated by upper and lower case [e.g., Q(I, J, K), q(i, j, k)].

1

(7)

3

2. Indices (i, j, k) correspond to the fine zone at the (left, bottom, back) of the coarse zone (I, J, K). 2

a)

b)

P(I,J,K)

c)

(I!1,J,K)

S 1(I,J,K)

#(i,j,k)

x3

Q

B3(I,J,K)

1

b3(i,j,k)

S1

"V

"V !V

x1

5

(I,J,K)

s 1(i,j,k)

"V

x2

B3

!V

!V

2

Figure 2: Restriction for: a) a volume-conserved zone-centred quantity; b) a volume-conserved 1face-centred quantity; and c) an area-conserved 3-face-centred quantity for a refinement ratio of 4.

Figure 1: The staggered mesh

S2 3 ! are restricted with B1 and B2 given by analogous expressions. Since B is area-conserved, “skin values” of B ! Note also that B ! satisfies the solenoidal condition numerically if !b does. (unlike volume-conserved S).

S2

The refluxing step

4

There are six cases, each indicated in Fig. 3. The first is for zone-centred volume-conserved quantities and is the same as standard AMR. The next three are for volume-conserved face-centred momenta and the last ! components, all peculiar to AMR on a staggered mesh. two are for area-conserved face-centred B 1. Q in a zone-centred RI immediately to the right of a fine grid is refluxed by: 



!−1 # q # !−1 # !−1 1  Q 1 ˜ f1 (i, j +β, k+γ, n+τ + 12 ) , F1 (I, J, K, N + 2 ) − Q(I, J, K) = Q(I, J, K) − ∆V β=0 γ=0 τ =0

(8)

where F1Q (I, J, K, N+ 21 ) is the coarse 1-flux of Q (units Q v ∆A ∆t) at the coarse 1-face (I, J, K) at coarse time step N + 12 (time-centred), and where f1q (i, j+β, k+γ, n+τ + 21 ) are the fine 1-fluxes of q at each of the %2 fine 1-faces contained within the coarse 1-face (I, J, K) at each of the % time-centred fine time steps.

&

1 F1S1 (I −1, J, K, N + 21 ) S˜1 (I, J, K) = S1 (I, J, K) − ∆V ) !−1 !−1 !−1 ' ( ### 1 ! ! s1 s1 1 1 f1 (i− 2 −1, j +β, k+γ, n+τ + 2 ) + f1 (i− 2 , j +β, k+γ, n+τ + 2 ) . − 2 γ=0 τ =0 β=0

The average between the fine fluxes at and arises by demanding that S1 integrated over the fine grid be the same as s1 integrated over the fine grid. 3. S2 in a 2-face-centred RI immediately and entirely to the right of a fine grid is refluxed by: 1 F1S2 (I, J, K, N + 12 ) S˜2 (I, J, K) = S2 (I, J, K) − ∆V −

!/2 !−1 !−1 # ##

β=−

where ηβ =

1 2

! γ=0 τ =0 2

)

ηβ f1S2 (i, j +β, k+γ, n+τ + 12 ) ,

(10)

for β = ± 2! , 1 otherwise.

4. S2 in a 2-face-centred RI immediately to the right and straddling the lower 2-boundary of a fine grid is refluxed by: 1 S˜2 (I, J, K) = S2 (I, J, K) − ∆V

&

% − 1 S2 F1 (I, J, K, N + 21 ) 2% −

!/2 !−1 !−1 # ##

ηβ f1S2 (i, j +β, k+γ, n+τ + 12 )

β=1 γ=0 τ =0

If the RI straddles the upper 2-boundary, the sum on β goes from − 2! to −1. 4

)

.



(12)

where E2 (I, J, K, N) is the coarse 2-emf (units B v ∆x ∆t) along the coarse 2-edge (I, J, K) at coarse time step N + 12 (time-centred), and where ε2 (i, j + β, k, n+ τ + 12 ) are the fine 2-emfs along each of the $ fine 2-edges contained within the coarse 2-edge (I, J, K) at each of the $ time-centred fine time steps.

5. 3-face-centred quantities require two 1-D interpolations within the 3-face: q(i+α, j +β, k) = Q(I, J, K) + δq1 (i+α) + δq2 (j +β).

i− !2

&

!−1 # !−1 #

˜3 (I, J, K) = B3 (I, J, K) + 1 E2 (I, J, K, N + 1 ) − B ε2 (i, j +β, k, n+τ + 21 ) , 2 ∆A β=0 τ =0

(9)

Figure 3: Refluxing replaces coarse fluxes/ emfs (black arrows) with suitable averages/ sums of fine fluxes/emfs (red arrows).

5. B3 in a 3-face-centred RI immediately to the right of a fine grid is “refluxed” by: 

2. S1 in a 1-face-centred RI straddling the right 1-boundary of a fine grid is refluxed by:

i− !2 −1

Figure 4: Ryu & Jones shock-tube problem 4a: (left) single fine grid solution; (middle) two fine grids (red) superposed on a coarse grid (black); and (right) % difference between the two solutions.

6. B2 in a 2-face-centred RI immediately to the right of a fine grid is “refluxed” by: 

!−1 # !−1 #



˜2 (I, J, K) = B2 (I, J, K) − 1 E3 (I, J, K, N + 1 ) − B ε3 (i, j, k+γ, n+τ + 12 ) . 2 ∆A γ=0 τ =0

(16)

This guarantees flux conservation in the 3-face. (13)

6. S3 is interpolated in the 3-direction by a linear interpolation between each pair q(i+α, j +β, k) and q(i+α, j +β, k+$). This ensures volume-conservation of S3 .

“Skin values” of B1 are restricted at the left and right boundaries of a fine grid, and are not refluxed there.

7. B3 is interpolated in the 3-direction, B2 in the 2-direction, and B1 in the 1-direction simultaneously in a manner suggested by Li & Li (2004, J. Comput. Phys., 199, 1) that enforces the solenoidal condition. Details are omitted for brevity.

Expressions for coarse zones to the left of a fine grid are obtained by adding 1 to the first index of F or E, adding $ to the first index of f or ε, and changing the sign in front of the factor 1/∆V or 1/∆A. Expressions for all other components at all other boundaries are obtained by a suitable permutation of indices and subscripts.

A test problem Fig. 4 shows the density profile of problem 4a from Ryu & Jones (1995, ApJ, 442, 228).

The reconstruction and restoration steps

The left panel shows the test done on a uniform grid with δx = 0.0025 (1,000 zones).

Interpolation of coarse grid variables to the fine grid (in constructing new grids or restoring boundary & = 0). values) is constrained by mass/momentum/flux conservation, and the solenoidal condition (∇ · B

The middle panel shows the same test on a computational domain with two fine grids (δx = 0.0025) fixed on a coarse grid (∆x = 0.005). The fast rarefaction is shown straddling a fine grid, the fast shock has passed through a fine grid, the slow shock has crossed the coarse-to-fine boundary and is now at the fine-to-coarse boundary, and the contact is at the coarse-to-fine boundary.

(11) A zone-centred quantity (e.g., ρ, eT ) is zone-centred in all three directions, while a 3-face-centred quantity (e.g., s3 , B3 ) is zone-centred in the 1- and 2-directions, staggered in the 3-direction. 1. To perform a 1-D interpolation in the 1-direction, start by computing the interface value, Q∗ (I), from a cubic fit to the data Q(I −2 : I +1). 2. Compute the parabolic interpolation function, Q1 (x1 ), for zone I that passes through Q∗ (I) and Q∗ (I +1) and satisfies the conservation constraint: Q(I, J, K)∆x1 (I) =

!−1 #

α=0

&

'

Q1 x1 (i+α) δx1

& ' # 1 !−1 ⇒ Q(I, J, K) = Q1 x1 (i+α) . r α=0

(14)

The interpolation parabolæ, Q1 (similar to those described by Colella & Woodward, 1984, J. Comput. Phys., 54, 174, but without the monotonisation and steepening steps), are continuous at zone faces. &

'

3. Set δq1 (i+α) = Q1 x1 (i+α) − Q(I, J, K). 4. Zone-centred quantities require three 1-D interpolation sweeps: q(i+α, j +β, k+γ) = Q(I, J, K) + δq1 (i+α) + δq2 (j +β) + δq3 (i+γ),

(15)

where α, β, γ = 0, 1, . . . $ − 1. This does not ensure monotonic interpolations and thus positivity. For positive-definite quantities, non-monotonic parabolic interpolations are replaced with linear (van Leer, 1977, J. Comput. Phys., 23, 276) interpolations which are monotonic even directionally split. 5

The right panel shows percentage errors (comparing non-AMR to AMR simulations). The largest errors (∼ 3%) are confined to zones within discontinuities which, by definition, are already in error by 50% or more from the analytical solution since shocks and contacts should have zero width. Elsewhere, errors are generally much less than 1%, indicating reflection and transmission errors of waves as they pass across grid boundaries are minimal.

A demonstration The background image shows a 2-D axisymmetric simulation of a proto-stellar jet computed from the subAU scale (bottom) where the jet is magneto-centrifugally launched from the proto-stellar accretion disc, to the 1,000 AU scale (top) where the jet is on the threshold of observable scales. It shows how the jet evolves from the knotty, magnetically-dominated regime at the smallest scales, to the hydrodynamicallydominated regime in which the jet exhibits the more familiar structures (e.g., bow shock, narrow jet) seen in observed proto-stellar jets (e.g., HH34). Shown are colour contours of the Alfv´enic Mach number (deep blue ⇒ MA = 0, aqua ⇒ MA = 1, black ⇒ MA > 100) with black magnetic field lines superimposed. The simulation shows the jet at t ∼ 100 yr. Jon Ramsey’s talk at 4:30 pm on Friday will give more details on this simulation.

Acknowledgements This project is supported financially by NSERC. Computational resources are provided by ACEnet which is funded by the CFI, ACOA, and the provinces of NS, NL, and NB. 6