Go to previous page Go up Go to next page

4.3 The magnetic-field divergence-free constraint

The numerical advantage of using the integral form of the conservation law, Equation (76View Equation), for the state vector variables does not apply in a straightforward way in the case of the GRMHD equations for the magnetic field components. While the induction equation is part of the hyperbolic system, there remains an additional constraint equation to be fulfilled, namely the magnetic field divergence-free constraint. Indeed, there is no guarantee that such divergence is conserved numerically when updating the magnetic field using the numerical procedure of Equation (76View Equation) for the rest of the components of the state vector (continuity, momentum, and energy equations).

Among the methods designed to preserve the divergence of the magnetic field (see [401Jump To The Next Citation Point] for a review), the constrained transport method (CT hereafter) designed by [119Jump To The Next Citation Point] and first extended to HRSC methods by [343] (see also [229] for a recent discussion), stands as the preferred choice for most groups. We note, however, that for AMR codes based on unstructured grids and multiple coordinate systems (such as those of [20Jump To The Next Citation Point286Jump To The Next Citation Point]) there are other schemes, such as projection methods or hyperbolic divergence cleaning, which appear more suitable than the CT method to enforce the magnetic field divergence-free constraint. The projection method involves solving an elliptic equation for a corrected magnetic field projected onto the subspace of zero divergence solutions by a linear operator. More precisely the magnetic field is decomposed into a curl and a gradient as

⃗∗ ⃗ ℬ = ∇ × A + ∇ φ , (95 )
whose divergence leads to an elliptic (Poisson) equation Δ φ = ∇ℬ⃗∗, which can be solved for the scalar quantity φ. The magnetic field is next corrected according to ⃗ℬ = ⃗ℬ∗ − ∇ φ. We note that an alternative projection scheme can be implemented by taking the curl of Equation (95View Equation) and solving for the vector potential A⃗.

Correspondingly, the basic idea of the hyperbolic divergence cleaning approach is to introduce an additional scalar field ψ, which is coupled to the magnetic field by a gradient term ∇ψ in the induction equation. The scalar field ψ is calculated by adding an additional constraint hyperbolic equation given by

d ψ 2 --- = − ch(∇ ⃗ℬ). (96 ) dt
Therefore, divergence errors are propagated off the grid in a wave-like manner with characteristic speed ch. Further details on these two methods are given, e.g., in [401Jump To The Next Citation Point].

The approach followed by the CT scheme to ensure the solenoidal condition of the magnetic field is based on the use of Stokes theorem after the integration of the induction equation on surfaces of constant t and i x, Σt,xi. Let us write the induction equation as

-1-∂-⃗ℬ ⃗ ⃗ √ γ-∂t = ∇ × Ω, (97 )
where we have defined the densitized magnetic field vector √ -- ℬ⃗ = γ⃗B and ⃗Ω = (α⃗v − ⃗β ) × B⃗.

Following [24Jump To The Next Citation Point] we can obtain a discretized version of Equation (97View Equation) as follows. At a given time, each numerical cell is bounded by six two-surfaces. Let us consider the two-surface Σ 3 t,x, defined by t = const. and 3 x = const., and the remaining two coordinates spanning the intervals from 1 x to 1 1 x + Δx, and from 2 x to 2 2 x + Δx. The magnetic flux through this two-surface is given by

∫ ⃗ ⃗ ΦΣt,x3 = Σ 3 B ⋅ dΣ. (98 ) t,x

Integrating Equation (97View Equation) on the two-surface Σt,x3, and applying Stokes theorem to the right hand side leads to the following equation

dΦ ∫ ---Σt,x3 = Ωidxi, (99 ) dt ∂(Σt,x3)
with ∂ (Σt,x3) being the boundary of the corresponding two-surface. This equation can be integrated to give
∫ t+ Δt∫ Φt+ΣΔt − ΦtΣ = ˆΩidxi dt, (100 ) t,x3 t,x3 t ∂(Σt,x3)
where the caret denotes that quantities ˆ Ωi are calculated at the edges of the numerical cells, where they can be discontinuous. At each edge, these quantities can be calculated using the solution of four Riemann problems between the corresponding faces, whose intersection defines the edge. It is important to emphasize, however, that irrespective of the expression used to compute ˆ Ωi, the method to advance the magnetic fluxes at the faces of the numerical cells satisfies, by construction, the divergence constraint (see, e.g., [24Jump To The Next Citation Point] for details). Indeed, the total magnetic flux through the closed surface S formed by the boundary of the numerical cells and enclosing a volume V can be calculated as
∮ ∫ Φ = ⃗ℬ ⋅ d⃗S = ⃗∇ ⋅ℬ⃗dV = 0 (101 ) T S=∂V V
after applying the Gauss theorem and divergence constraint. This equation implies that no magnetic flux source exists inside the volume V and, therefore, the magnetic flux is a conserved quantity as dΦT ∕dt = 0. As a result, every numerical scheme constructed on the basis of the CT scheme will conserve magnetic flux up to machine accuracy. If one is able to generate initial conditions, which satisfy the divergence constraint, i.e., with Φ = 0 across the boundary of each numerical cell, then such a constraint will be preserved during the numerical evolution.

Any of the Riemann solvers and flux formulae discussed in Section 4.1.2 can be thus used to calculate the quantities Ωˆi needed to advance in time the magnetic fluxes following Equation (100View Equation). At each edge of the numerical cell, ˆΩi is written as an average of the numerical fluxes calculated at the interfaces between the faces whose intersection define the edge. Let us consider, for illustrative purposes, ˆ Ωx. If the indices (j,k,l) denote the center of a numerical cell, an x −edge is defined by the indices (j,k + 1∕2,l + 1 ∕2). By definition, y z z y Ωx = α (&tidle;v B − &tidle;v B ). Since F y(Bz ) = v&tidle;yBz − &tidle;vzBy and F z(By ) = &tidle;vzBy − &tidle;vyBz, we can express ˆΩx in terms of the fluxes as follows

ˆ 1-ˆy ˆy Ωx j,k+1∕2,l+1∕2 = 4[Fj,k+1∕2,l + Fj,k+1∕2,l+1 ˆz ˆz − F j,k,l+1∕2 − F j,k+1,l+1∕2], (102 )
where ˆF y(ˆFz) refers to the numerical flux in the y (z) direction corresponding to the equation for z B (y B) and multiplied by α to account for the correct definition of Ω. Analogous expressions for Ωˆy j+1∕2,k,l+1∕2 and ˆΩz j+1∕2,k+1∕2,l, can be derived in a similar way. We refer the interested reader to [401] for additional details on the constrained transport method.
  Go to previous page Go up Go to next page