FIG. 1: Demonstration of the effectiveness of divergence cleaning at handling deviations from a divergenceless mag- netic field. A calculation of the divergence of the magnetic field, x 2 [ ⃗ ∇ · ⃗ B ] , along the x -axis for three times is shown for each of three cases: (top) No divergence cleaning; (middle) Cleaning with c h = 0 . 1 and c r = 0 . 01; and (bottom) Cleaning with c h = 1 . 0 and c r = 0 . 1. These runs were otherwise iden- tical for a magnetized, rotating star of coordinate equatorial radius of 10 with spin along the z -axis with three levels of re- finement (at ± 50, ± 25, and ± 12 . 5) and a perturbation to the magnetic field to introduce an explicit deviation from diver- genceless. Concentrating on the MHD equations, the metric terms were frozen at their initial values, e.g. the Cowling ap- proximation. The bottom row shows clear wavelike behavior as it “cleans” the divergence. can be solved for P by substituting in ρ$_{0}$ = D/W and the EOS to obtain Combining these equations, f is a function of a single unknown x . This equation can be solved numerically using the Newton-Raphson method. It is useful to note that x has a minimum physical value, which is found by requiring in Eq. (33) that W 2 ≥ 1. The code is constructed within the had computa- tional infrastructure which provides distributed adaptive mesh refinement (AMR). The AMR follows in the style For initial data generated with Magstar , the diver- gence is around machine precision, and so to test the implementation of divergence cleaning, we introduce a perturbation to the magnetic field in order to produce a significant amount of divergence to the magnetic field. In The Maxwell equations require that the magnetic field be divergenceless. This is the so-called “no monopole” constraint. A variety of schemes exist with the goal of controlling the growth of the divergence. We choose a strategy that ensures flexibility and robustness when dealing with multiple grid structure (as in AMR) and al- lows, in principle, for a clean boundary treatment [44]. To this end we have implemented hyperbolic divergence cleaning as described in [36] (also see [45]). We thus introduce a scalar Ψ( x, y, z, t ) which is sourced by the negative of the divergence of the magnetic field as shown in Eq. (30). As described in [36], this scheme implies that the divergence obeys a damped wave equation so that constraint violations propagate off the grid and their value is reduced. We have, however, found certain issues with precision occurring within the primitive solver. Typical maximum values of the density are about 10 − 2 in geometric units, with a floor value of δ = 10 − $^{8}$. We have found it use- ful therefore to scale Newton’s constant G such that the fluid densities and pressures are close to order unity. Thus, rather than using the typical choice of G = 1, we might use G = 1 / 1000. As G affects only the coupling of the fluid to the geometry, the evolution of the geomet- ric equations is not affected by this scaling. Empirically, we find that scaling G allows the primitive variables to be more easily recovered in low density regions. This im- provement appears to be related to finite precision effects in the primitive solver. The scaling decreases the effec- tive floor value while avoiding the problems associated with having a true vacuum. where δ is chosen to be many orders of magnitude smaller than the maximum densities and pressures in the initial data. The comparison of otherwise identical runs but with different floor values suggest that the use of an at- mosphere generally does not affect accuracy. The fluid equations diverge as the density goes to zero, and, as is standard practice, we disallow any true vacuum by setting such regions to a floor or atmosphere value. The floor is applied after each fluid update as of Berger & Oliger [38], but uses the tapering condi- tion for AMR boundaries instead of temporal interpo- lation [39]. The combined set of geometric equations and fluid equations, Eqs. (5-7) and Eqs. (26-30) respectively, is discretized using the method of lines. The geometric equations are discretized using operators that satisfy a summation by parts property [40, 41]. The fluid equa- tions are discretized using the HLLE method [42]. The semi-discrete equations are solved using a third order ac- curate, total variation diminishing (TVD) Runge-Kutta solver [43]. 5 D ← min( D, δ ) , τ ← min( τ, δ ) , (37) P = Γ − 1 Γ ( x W 2 − D W ) . (36) III. IMPLEMENTATION DETAILS