{'fraction_non_alphanumeric': 0.028593233273976087, 'fraction_numerical': 0.011727295853472399, 'mean_word_length': 4.661146313364055, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Axisymmetric Hydrodynamics in Numerical Relativity Using a Multipatch Method

We describe a method of implementing the axisymmetric evolution of general-relativistic hydrodynamics and magnetohydrodynamics through modification of a multipatch grid scheme. In order to ease the computational requirements required to evolve the post-merger phase of systems involving binary compact massive objects in numerical relativity, it is often beneficial to take advantage of these system's tendency to rapidly settle into states that are nearly axisymmetric, allowing for 2D evolution of secular timescales. We implement this scheme in the Spectral Einstein Code (SpEC) and show the results of application of this method to four test systems including viscosity, magnetic fields, and neutrino radiation transport. Our results show that this method can be used to quickly allow already existing 3D infrastructure that makes use of local coordinate system transformations to be made to run in axisymmetric 2D with comparable results.


I. INTRODUCTION
The detection of the gravitational wave signal resulting from the merger binary black hole systems by the LIGO and VIRGO collaborations [1][2][3][4][5][6][7] along with the detection of simultaneous electromagnetic and gravitational wave signals from binary neutron star mergers [8][9][10][11], and the corresponding need for theoretical predictions with which to compare them, has given renewed urgency to the goal of accurately modeling these systems throughout the merger process. For systems involving at least one neutron star, it is the post-merger state that is primarily responsible for the observable electromagnetic signals. Modeling of these systems through numerical relativity simulations provides critical insight into the dependencies of the signals on binary parameters and nuclear physics. Unfortunately, running these simulations in the high-resolution required to get accurate predictions can present large computational resource barriers in simulated time or size scales. However, the post-merger environment has a useful property: by taking advantage of these systems' tendency to approach an axisymmetric state, we can ease the computational resources required to simulate these systems over extended scales of both time and space. Although the dynamical timescales of remnant neutron stars and accretion disks, of the order ∼ms at most, are reasonably accessible to 3D simulations, secular effects that drive the subsequent evolution can operate on much longer timescales. Particularly important are angular momentum transport effects that can act on a wide range of timescales of up to hundreds of milliseconds [12,13], and neutrino cooling effects that operate on timescales of up to several seconds [14].
The use of axisymmetry in numerical relativity simulations has been explored by several groups. This typically involves evolving Einstein's equations using the cartoon method [15] while evolving hydrodynamics by writing the relevant equations in a cylindrical coordinate system [16][17][18]. The cartoon method does involve some loss of accuracy due to interpolations required in the method, and some effort has been made to avoid these [19]. Additionally, evolution problems due to the coordinate singularities that arise from the use of polar coordinate systems have been avoided by use of a reference metric [20][21][22]. Methods also exist that help with issues of spatial resolution on large scales, such as adaptive mesh refinement [19,23], which is able to concentrate resolution where it is most needed, while in most cases still building the grid from Cartesian domains. In multipatch methods [24][25][26][27][28][29], one introduces coordinate patches, each with its own local coordinate system in which it takes a simple shape (e.g. a Cartesian block), but which can be deformed in the global coordinate system and fit together into a grid to match the geometry of the problem. A number of methods used in numerical relativity not usually called "multipatch" have local coordinate systems and therefore fit into this general category, including the multidomain pseudospectral sector of the Spectral Einstein Code [30] and the multielement discontinuous Galerkin methods [31,32] which many hope will form the basis of the next generation of numerical relativity codes.
In this paper, we describe a method of implementing the axisymmetric evolution of the general-relativistic equations of ideal radiation hydrodynamics and magnetohydrodynamics through modification of a multipatch grid scheme, applicable to any method using the local patch coordinates framework, which we implement in the Spectral Einstein Code (SpEC) [30]. We show that this method is able to use the well-developed methods used for full 3D simulations to quickly be made to run in 2D axisymmetry.
In Sec. II, we describe the evolution equations for our (magneto)hydrodynamic variables and the application of our axisymmetry method to them. In Sec. III several tests of this axisymmetry method are presented: a stationary TOV star, a viscous differentially rotating star, a magnetized accretion disk, and neutrino radiation in a spherically symmetric supernova collapse profile, each showing good agreement with 3D results or previous axisymmetric simulations. Concluding remarks are given in Sec. IV, where we summarize our results and discuss future plans.

A. Evolution Equations
We use SpEC to evolve Einstein's equations and the general relativistic equations of ideal radiation (magneto)hydrodynamics. SpEC evolves Einstein's equations and the general relativistic hydrodynamics equations on two separate computational grids. A multidomain grid of colocation points is used to evolve Einstein's equations pseudospectrally in a generalized harmonic formulation [33] while the general relativistic (magneto)hydrodynamics equations in conservative form are evolved on a finite difference grid. The finite difference grid uses an HLL approximate Riemann solver [34]. Reconstruction of values at cell faces from their cellaverage values is done using a high-order shock capturing method, a fifth-order WENO scheme [35,36]. Time evolution is performed using a third-order Runge-Kutta algorithm with an adaptive time-stepper. At the end of each time step any necessary source term information is then communicated between the two grids, using a thirdorder accurate spatial interpolation scheme [37].
The following sections make use of the 3+1 decomposition of the spacetime metric where α is the lapse, β i the shift, and γ ij is the threemetric on a spacelike hypersurface of constant coordinate t. The three-metric is the projection onto spatial hypersurfaces of the four-metric: where n µ = (−α, 0, 0, 0) is the unit normal to the t = constant hypersurface. Additionally, we use the units with G = c = 1 throughout.

Fluid
We begin by treating our fluid as a perfect fluid with the stress-energy tensor where ρ 0 is the baryon density, h = 1 + P/ρ 0 + is the specific enthalpy, P is the pressure, u µ the four-velocity, and the specific internal energy. The general relativistic hydrodynamics equations are evolved using the conservative variables where W = 1 + γ ij u i u j is the Lorentz factor and γ is the determinant of γ ij . Using conservation of energy and momentum, ∇ ν T µν = 0, and baryon number conservation, ∇ µ (ρ 0 u µ ) = 0, we get the evolution equations for the conservative variables: where the Eulerian velocity v i is related to the fluid transport velocity v T i by v T i = αv i − β i . Additionally, for simulations involving nuclear matter and neutrinos, we evolve the electron fraction of the fluid, Y e , To close these equations we must also supply an equation of state for the pressure and enthalpy: P = P (ρ * , T, Y e ) and h = h(ρ * , T, Y e ).

Magnetic Fields
To handle magnetic fields, we begin by adding the electromagnetic contribution, T EM µν , to the fluid stressenergy tensor, where and F µν is the Faraday tensor. We treat the fluid as a perfect conductor, F µν u ν = 0, which gives a fixed electric field for a given magnetic field. We use two different methods for evolving the magnetic field, as described in [29]. The first method evolves the magnetic vector potential A i and scalar potential Φ. In the generalized Lorentz gauge [38], the most robust gauge choice we have explored, the evolution equations are where ξ is a specifiable constant of the order of the mass of the system. The second method evolves the magnetic field using a covariant hyperbolic divergence cleaning method [39][40][41] in which an auxiliary scalar evolution variable Ψ is introduced in order to propagate and damp monopole formation. In this method, the induction equation takes the form whereB i = √ γB i ,Ψ = √ γΨ, K i i is the trace of the extrinsic curvature, and λ is a specifiable damping constant.

Neutrinos
Neutrino evolution is handled using the gray twomoment scheme as described in [42,43]. This method provides evolution of neutrino average energy densities, flux densities, and number densities. We define three neutrino species that we evolve: electron neutrinos ν e , electron antineutrinosν e , and the heavy lepton neutrinos ν x . The heavy lepton neutrino species groups together the four heavy lepton neutrinos and antineutrinos: ν µ , ν µ , ν τ , andν τ .
We can describe each of our three species of neutrinos ν i using each species' distribution function f ν (x µ , p µ ), where x µ = t, x i gives the time and position of the neutrinos and p µ is the 4-momentum of the neutrinos. f ν evolves in phase space according to the Boltzmann transport equation: where the term C f (ν) includes all collisional processes (emissions, absorptions, and scatterings). We simplify the radiation evolution by taking the gray approximation (integrating over the neutrino spectrum) and evolving the lowest two moments of the distribution functions of each neutrino species, truncating the moment expansion by imposing the Minerbo closure [44]. Our evolved quantities are projections of the stressenergy tensor of the neutrino radiation, T rad µν . The decomposition of T rad µν in the fluid frame is with H µ u µ = S µν u µ = 0. The energy density J, flux density H µ , and stress density S µν of the neutrino radiation as observed in the frame comoving with the fluid are related to the distribution functions by where ν is the neutrino energy in the fluid frame, dΩ denotes integrals over solid angle in momentum space, and where l α u α = 0 and l α l α = 1. We also make use of the decomposition of the neutrino radiation stress-energy tensor as observed by a normal observer, with F µ n µ = P µν n µ = F t = P tν = 0. Additionally, for each species of neutrino we consider the number current density: where N is the neutrino number density, and F µ is the number density flux. The decomposition of N µ relative to the fluid frame can be expressed in terms of J, H µ , and the fluid-frame average neutrino energy ν as We define a projection operator onto the reference frame of an observer comoving with the fluid, This allows us to then use the fluid-frame variables to write equations for the energy, flux, and stress tensor in the normal frame (i.e. the frame with 4-velocity equal to the normal vector) by making use of the decomposition of the 4-velocity, can then be written in conservative form: whereP ij = √ γP ij . Complete treatment of these equations requires prescriptions for the closure relation that computes P ij (E, F i ), the computation of F j and the collisional source termsS α rad and C (0) which couple the neutrinos to the fluid (and introduce corresponding source terms to the right hand side of Eq. (8), (9), and (10)). Details on the treatment for these are beyond the scope of this paper, and are available in [42,43].

Viscosity
Viscosity is implemented using the approach of [45] that extends the Newtonian large-eddy simulation framework to general relativistic systems. In the large-eddy simulation framework, we recognize that although the equations for energy and momentum evolution allow for evolving modes at all scales, in numerical simulations on a discrete grid we can only evolve modes for which we have sufficient resolution to cover. Thus each computational cell deals with averaged values, while any modes smaller than the cell are removed.
We therefore average over and filter out small scales in the velocity field, leaving equations for the resolved fields: where K ij is the extrinsic curvature and S ij = S i v j + P γ ij . In order to complete this set of mean-field equations, we must provide a closure condition for the quantity S i v j : τ ij is the subgrid scale stress tensor, that captures the turbulent modes unresolved by our grid. We model this tensor using where ∇ is the covariant derivative compatible with γ ij . The quantity ν T possesses a dimension of a viscosity, which leads us to the assumption where c s is the sound speed of the local fluid. mix is the characteristic length over which our subgrid scale turbulence occurs and is known as the mixing length.
As explained in [46], we find that, to maintain the relations Eq. 4-6 for resolved fields, Eq. 32 must be altered. In this paper, we use the energy equation with the correction to 2nd order in v, which is

B. Multipatch Axisymmetry
Multipatch methods work by dividing the computational domain into separate domain patches, each of which may have its own local coordinate system x i L related to the global coordinate system x i G by a map which controls the embedding of the domain in global space. In local coordinates, the patch is (for all applications in this paper) a simple Cartesian grid. The basis vectors ∂/∂x i L and ∂/∂x i G are then related by the Jacobian transformation matrix of the map. Importantly, the patches may have differing shapes in the global coordinate system that can be tailored to better capture the desired features of the simulation. Since our evolution equations for the conservative variables are generally covariant, evolution can be performed directly in the local coordinate system of each individual patch and then the result can be transformed back to the global coordinate system for any necessary communication of information between patches.
Communication between domain patches occurs through synchronizing values in the ghost zones of each patch at the end of each timestep. In the case that these subdomain patches overlap but do not have directly matching points we communicate data by interpolating values between points. Additionally, we create ghost zone points that extend beyond any symmetry boundaries that we have defined in order to impose boundary conditions. During the communication phase, these ghost zone points are filled with data from the live points using the appropriate symmetry conditions (i.e. axisymmetry or a reflection symmetry).
When evolving a three-dimensional system using a two-dimensional computational domain, each gridpoint represents a ring labeled by two nonazimuthal coordinates. Quite general 2D maps are possible to relate local to global coordinates, but two are particularly useful. A linear map (x i G = a i x i L + b i ) corresponds to patches that are globally rectangular blocks, covering cylinders in 3D. A polar map [e.g.
] corresponds to patches that are globally wedges of circles, covering a specified range of polar r, θ. A combination of wedges covering 0 < θ < π in 2D covers a spherical shell domain in 3D. A general 2D grid can contain arbitrary combinations of rectangular blocks and wedges, as shown in Fig. 1.
Although the grid is 2D, the tangent space on which vectors live is still 3D; even axisymmetric systems can have azimuthal velocity and magnetic field components, for example. The third coordinate in the local coordinate system is set to be the global azimuthal φ. Then the local coordinates for a rectangular block will be (up to linear transformation) cylindrical-polar, while the local coordinates for a wedge patch will be (up to linear transformation) spherical-polar. By modifying the map Jacobian, we can make the existing transformation between local and global coordinates handle transforming the third coordinate into an azimuthal coordinate that can be used to perform axisymmetric evolutions. To do this we expand the elements of the Jacobian matrix using the chain rule to add in the effects of the polar transfor-mation: where x G are the global coordinates, x L are the local coordinates of a given grid patch, and x A are a set of global polar coordinates. Since the global and polar coordinates only differ in terms involving the azimuthal direction, the final change from the original Jacobian, J i j , to the new axisymmetry Jacobian, J axi i j , will be straightforward: where is the coordinate distance from the rotational symmetry axis and we have chosen coordinate directions 1 and 2 to correspond to the two coordinates defined by our two-dimensional computational domain and coordinate direction 3 is transformed to the axisymmetric azimuthal direction φ. We also make use of the Hessian matrix in the transformation of the derivatives of metricrelated quantities to the local coordinates, and must likewise make similar adjustments to the Hessian: Explicitly, Generally the evolution of Einstein's equations using SpEC's pseudospectral grid tends to use much less computing time than the hydrodynamics evolution, so our axisymmetry method is primarily aimed at implementing axisymmetric evolution on the hydrodynamics grid while evolving Einstein's equations in 3D. Information required by the pseudospectral grid from the hydrodynamics grid is expanded back to 3D during communication. We mention that, for spherical shell pseudospectral domains, whose colocation points correspond to an expansion of functions in terms of spherical harmonics, azimuthal information can be reduced by reducing azimuthal resolution, corresponding to a lowering of the azimuthal mode number m retained in spectral expansions. It cannot be lowered to m max = 0 because the spectral evolution uses Cartesian components of tensors. We find, however, that the speed increase from doing so is modest, and the resulting spectral grids are more prone to constraint-violating instabilities, so we have not used azimuthal resolution reduction on pseudospectral grids for the simulations in this paper.
The conservative form of radiation magnetohydrodynamics evolves variables that are densities and thus proportional to √ γ. Under local to global transformation, the metric determinant transforms as where J is the determinant of the Jacobian. Note that J is zero on the axis, and indeed would naturally change sign there because the orientation of the basis vectors switches there. SpEC always takes a positive square root, but the only points on the other side of the axis are ghost zone points (needed to impose the symmetry boundary conditions), and non-smooth functions like √ γ are not interpolated or reconstructed. Unfortunately, when evolving, this method is prone to producing errors near the symmetry axis that, without correction, grow over time. Vector and tensor valued quantities are most heavily affected due to direct transformation of components in the azimuthal coordinate direction introducing singular terms. An example of this type of error is shown in Fig. 2. Eventually though, all of our evolved quantities, including scalar quantities, will suffer from errors due to also picking up a singular term in the determinant of the 3-metric.
The problem primarily occurs during the computation of the divergence of the flux term, F A , in the evolution equation of a given quantity A with S A being any source terms appearing on the righthand side of the equation. Some early 2D general relativistic hydrodynamic simulations stabilized the axis evolution using dissipation [47,48]. Our solution, inspired by [22], is to factor out singular terms that have been introduced to F A during the transformation to the local coordinates prior to computing the divergence. Depending on the specific component of the flux F A corresponding to A, there may be multiple factors of that need to be removed: whereF A is just the -factored form of the flux, and the integer n will depend on A. We can now instead take the divergence of this factored form of the flux and apply the chain rule, which gives We can also take advantage of the property that if the coordinate specified by corresponds to one of the directions in the global coordinate system, for example if the global coordinates are Cartesian, the derivatives of with respect to the local coordinates can be directly taken from components of the Jacobian dealing with the direction associated with . With this, all of the singular terms introduced from the polar Jacobian are removed from the divergence. Importantly though, the divergence ofF A in the first term on the right side of this equation will need to be computed using the value ofF A at cell faces using the Riemann solver, whileF A in the second term on the right side will use the value at cell centers.
Additionally, since all components have now been transformed into a polar coordinate system, from the definition of axisymmetry we have where the φ-index indicates the coordinate of the axisymmetric azimuthal direction. This allows us to ignore the azimuthal portion of the divergences so that we only need to apply the factoring to the two components of the flux that lie in the plane of the computational grid (i = 1 and 2 in the below factoring). All of our evolved quantities carry a factor of √ γ, which will also acquire a singular term, from the transformation of γ ij to the local coordinate system, that also needs to be handled analytically. The flux factoring thus falls into three broad categories for our current evolution equations. Factoring of fluxes for scalar density quantities [Eq. (7), (8), (10), (13), (15), (29), (31), and the added term in (37)], takes the form Factoring for covariant vector density quantities [Eq. (9), (12), and (30)], takes the form Factoring for contravariant vector density quantities [Eq. (14)], takes the form In each of these, the index i only covers coordinates 1 and 2 due to Eq. 47. SpEC and most other relativistic hydrodynamics codes use conservative shock capturing techniques with approximate Riemann solvers. For codes of this type, a convenient way to implement this factoring program is to use a different coordinate basis, with 1 ∂ ∂x φ instead of ∂ ∂x φ , on cell faces than on cell centers. That is, one simply reconstructs factored quantities. When evolving a magnetic vector potential, it is also necessary to factor A φ when computing B i .
Metric-related quantities (γ ij , α, β) are evolved on their own separate spectral grid in 3D and are communicated to the hydrodynamics grid at the end of each time step. Spatial derivatives of these metric quantities are computed while on the metric grid and then communicated to the hydrodynamics grid, at which point they can be transformed into the local coordinate system as needed. The transformation to local coordinates uses the analytic Jacobian and Hessian, so metric derivatives automatically have their singular factors treated analytically. The transformation equations for global to local components of metric derivatives are

C. Auxiliary Entropy Variable
After each substep, the evolved variables (ρ * ,τ ,S i ,ρ * Y e ,B i ) must be used to recover the primitive variables (ρ 0 ,T ,Y e ,u i ,B i ), a process that involves multi-dimensional root-finding. In particular, if the internal energy is small compared to kinetic or magnetic energy, the temperature recovered from total energy and momentum densities will be unreliable. Due to numerical error, recovered T and u i , especially at very low densities, may be unphysical, or there may not even be a set of primitive variables corresponding to the evolved variables at a point.
As in [29,50], we introduce an auxiliary entropy density evolution variable ρ * S, where S is the specific entropy. The variable ρ * S obeys a continuity equation (viscous and neutrino source terms being unimportant for its purpose) which can be treated in axisymmetry like the other scalar density evolution equations. After each substep in time, SpEC first attempts to recover primitive variables using the standard evolution variables. If this is not possible, or if the recovered specific entropy decreases by more than a fixed percentage compared to its advected value [51], primitive variables are recovered disregarding τ and using ρ * S. At the end of each substep, the primitive variables are used to reset all evolution variables, so that τ and ρ * S are synchronized to each other.
For physical equations of state (e.g. finite-temperature nuclear-theory based EoS), the actual statistical mechanical entropy per baryon can be used to define S. However, in numerical relativity, equations of state are commonly used which have no uniquely defined entropy or temperature, although with absolute zero specified from outside (e.g. for Gamma-law EoS, a value of the polytropic constant is defined to be "cold"). A common case is an EoS with nuclear physics-motivated cold component plus a simple thermal Gamma-law component added on.
In terms of baryonic number density n = ρ 0 /m amu and internal energy density u, where The first law gives Combining the three above equations yields, after a short calculation, so (u−u c )ρ 0 −Γ th advects for adiabatic change, indicating that this is an acceptable S variable. For Gamma-law EoS, one can set u c = 0, yielding the standard auxiliary entropy variable (up to a scaling factor) for this case. . We see an initial spike in density at the center of the star. This is caused by relaxation of the surface of the star, creating a disturbance that moves inward. Density is continuous but not smooth at the surface, so this featuret exhibits approximately first-order convergence. After this initial peak settles, we see second-order convergence.

B. Differentially Rotating Star
We choose a star with very similar profiles and global properties as the differentially rotating star used in [52] and likewise use this star to test the evolution of a system under the influence of viscosity. The star has an initial baryonic rest mass of 2.64 M and an equatorial radius R e = 10.2 km. We use a piecewise polytropic equation of state, in two pieces, of the form where κ 1 and κ 2 are polytropic constants, Γ 1 and Γ 1 are the polytropic indices, and ρ t is the density at which we transition between the two pieces. For this star, we choose the polytropic indices to be Γ 1 = 4/3 and Γ 2 = 11/4; we set the transition density between the two to be ρ t = 1.91 × 10 14 g cm −3 and the low-density polytropic constant to κ 1 = 0.15GM 2/3 . The initial rotation profile for the star is given by u t u φ =Â(Ω 0 − Ω) where Ω 0 is the angular velocity along the rotation axis and we chooseÂ = 0.8R e . The ini- tial equilibrium state is supplied by the code of Cook, Shapiro, and Teukolsky [53].
In order to handle outflows that will occur when viscosity is enabled, we create a computational grid better suited for resolving both the central star and low density outflowing material. Since any outflows that occur will rapidly drop in density and are not expected to have any small detail features of concern after they leave the region of the star, we leverage the utility of the multipatch technique to apply differing grid structures to each zone of interest. In the central region containing the star we employ the same rectangular grid structure as seen in the previous TOV star test, with a resolution of 100×100 grid points. In the outflow region we switch to a polar grid with constant latitude resolution (so that the proper spacing between angularly adjacent points increases with distance from the star). The polar grid has 50 points in the angular direction (covering 0 < θ < π/2) and 400 points in the radial direction. We apply a map to the entire grid that allows us to reduce radial resolution at large distances: where r is the radius in grid coordinates (the coordinates in which radial grid spacing is uniform), and R is the radius in the original quasi-isotropic, asymptotically-Minkowski coordinates. The map provides an approximately linear grid spacing for radii less than β, which we have chosen to be at 25.85 km, and then switches to an exponential grid spacing based on γ, which is chosen such that r outer = 73.5 km is mapped to R outer = 2205 km. The pseudospectral grid used for the evolution of Einstein's equations is composed of an inner ball at the center of the star surrounded by a series of spherical shells extending to a distance of 2940 km. We impose a density floor outside of the star which is necessary to avoid division by zero in our finite difference solver. At densities below the floor we recover temperature and velocity using the prescription described in [54]. We recover the primitive variables from the conservative variables using the auxiliary entropy variable in these areas using the process in [29]. In this test, we have modified the density floor from our previous implementations to use a floor dependent on radius: where we have chosen A = 1.62 × 10 4 g cm −3 and B = 1.62 × 10 −2 g cm −3 . For this test, we employ the viscosity treatment described in Section II. To make a comparison with the results of the α-viscosity model used in [12], we devise a mixing length mix corresponding to the same kinematic viscosity as a constant α. The α-viscosity model is generalized to differentially rotating stars in [12] by setting where c s is again the local sound speed and Ω e is the angular velocity of the star at the surface on the equator. Equating this to Eq. (36), we can get an approximate relation between the strength of a given mixing length to that of an α-viscosity parameter: As evolution begins the star quickly begins to transport angular momentum outward causing the rotational velocity profile to become flatter. Although the rotation profile does flatten, we see from Fig. 5 that the profile never completely settles into a rigidly rotating state, and retains some differential rotation. This is a feature of this viscosity method [46].
Qualitatively the outflow near the star produces the expected distribution of material producing a short, low density burst of material as viscosity is enabled, and at later times as more material leaves the star a disk begins to form.

C. Magnetized Disk
We evolve a standard axisymmetric MHD test problem: a magnetized torus around a Kerr black hole. The initial conditions for this test are matched to the "fiducial model" of McKinney and Gammie [55]. A black hole with dimensionless spin J/M 2 = 0.938 is surrounded by a Fishbone-Moncrief torus [56] with inner edge at r BL = 6M and specific angular momentum determined by u t u φ = 4.281. The torus has initial maximum density ρ 0 = 1 and a Γ = 4/3 equation of state. A confined poloidal seed field is introduced via the initial vector potential 1-form A = A 0 max(ρ 0 − 1, 0) dφ, with A 0 chosen to make the maximum ratio of magnetic to gas pressure be around 0.01. We evolve for 3000M on a 256×256 spherical-polar grid with inner radius at r BL = 1.32M and maximum radius at 60M .
The Kerr spacetime is written in in Kerr-Schild coordinates. We make the standard change of variables for spherical-polar disk simulations: Setting a uniform grid in x 1 , x 2 concentrates resolution near the black hole and on the equator. We set h = 0.5. Finally, because r = r BL we compose with a final coordinate map to map the coordinate spheres (x 2 +z 2 ) 1/2 = C to surfaces of constant Kerr radius r BL = C. This allows an excision inner boundary inside the horizon r BL = r + that conforms better to the horizon shape.
For this run, we use a position-dependent density floor ρ 0 > 10 −5 r −3/2 . We also increase ρ 0 and P in the magnetically-dominated region as needed to maintain b 2 /ρ 0 < 10 and b 2 /P < 500, which significantly improves the step size chosen by the adaptive timestepper. We evolve both with hyperbolic divergence cleaning and vector potential evolution. For the vector potential evolution, we use the generalized Lorentz gauge [38]. Simpler gauges, such as the algebraic ∂ t A = v · B and advective ∂ t A = −L v A give the same evolution of gauge-invariant quantities but, after a while, at a drastically reduced timestep, presumably because the vector potential does not remain as smooth.
The vector potential evolution benefits from added explicit dissipation. We apply Kreiss-Oliger dissipation [57] to the evolution of A i and Φ with a coefficient of 0.001. (Our dissipation operator is defined as a sum of fourth derivatives with respect to local coordinates but applied to global components of the relevant evolved variables.) Without dissipation, grid-scale ripples appear in the magnetic field atop an otherwise reasonable field structure. If the coefficient is increased to 10 −2 , the main difference is a slightly lower asymptotic speed in the polar jets. Kreiss-Oliger dissipation is not needed for divergence cleaning runs; in fact, it destabilizes the magnetic field evolution near the excision zone. Instead, extra dis- sipation for divergence cleaning simulations is obtained by setting the maximum signal speeds in our HLL Riemann solver for the evolution ofB and Ψ to the null speeds.
The qualitative expectations for this problem are wellknown and are reproduced for our runs for both types of B field evolution. Magnetic winding generates a toroidal magnetic field, while the magnetorotational instability triggers turbulence in the disk. Matter falls into the black hole at an average rate of aboutṀ ≈ 10 −1 . The poles become magnetically dominated. An outgoing Poynting flux can be found in this region, and gas accelerates to near the speed of light on the poles away from the black hole. The magnetic field energy grows for the first 1000M , then saturates, then begins to die away at a steady rate. This decrease of the magnetic field is not physical but it is expected in any axisymmetric simulation (at least one not enhanced by dynamo-modeling additions to the induction equation [58]) because of the anti-dynamo theorem. Outside the region close to the poles, a mildly relativistic wind is seen. The configuration of the system at t = 1500M is shown in figure 6.
None of this is newsworthy, although it is reassuring to confirm for the first time that SpEC can produce magnetically-dominated jets when they are expected. For our purposes, the main value of this test is that we can check, for a complex, astrophysically interesting MHD problem, that our code produces no unphysical axis artifacts in any quantity we have checked (ρ 0 , v i , B i , b 2 /P ). Of course, the axis actually is a special region in this problem, which is clearly seen in the solution, but this can easily be distinguished from artifacts of the coordinate singularity because the latter have gridspacing width. The absence of such glitches is, in fact, a nontrivial accomplishment. For divergence cleaning evolutions without factoring of the evolution equations, grid-scale axis artifacts in the velocity are easily seen, although they can be suppressed by using low-order reconstruction (MC2 [59]) near the axis. For vector potential evolutions without factoring the computation ofB from A, axis glitches become so severe that simulations crash shortly after accretion onto the black hole begins.
Although the results are qualitatively similar, we consider the vector potential method superior for this problem, at least with our current implementations. In divergence cleaning methods, Ψ builds up at boundaries, particularly the excision boundary. The amount tends to grow with time and we fear would eventually endanger the simulation. Because of it, magnetic energy fluxes are not reliable in the inner layer of points (while in vector potential evolutions, the inner layer shows no problems). Presumably the solution would be to improve the treatment of the magnetic variables at boundaries.

Spherically Symmetric Collapse Profile
Initial testing of the axisymmetric neutrino code was performed by comparing the results obtained from the spherically symmetric post-bounce supernova profile used in [43] in both 2D axisymmetry with equatorial symmetry and in 3D using octant symmetry. In this test we evolve the moments of the neutrino distribution function, fluid temperature, and fluid composition (the electron fraction Y e ) for a 1D profile constructed as a spherical  ν e 2D ν a 2D ν x 2D ν e 3D ν a 3D ν x 3D FIG. 8. Comparison of neutrino number density along the z-axis for the 2D and 3D spherically symmetric supernova collapse profile at t = 1.5 ms. νe, νa, and νx represent the electron neutrinos, electron antineutrinos, and heavy lepton neutrinos respectively. average of a 2D core collapse simulation 160 ms after bounce. The velocity of the fluid is set to zero.
We perform this test in 2D on a square grid with length 300 km and a resolution of 200 × 200 grid points. In 3D, we use a cube with the same length of 300 km and a resolution of 200 × 200 × 200 grid points. Both systems were evolved for 1.5 ms using a fixed timestep of 4.9 × 10 −3 ms to ensure that no error was introduced from possible differences between 2D and 3D in the adaptive timestepper. The 2D test used 24 processing cores and required 2.41 core-hours of run time, whereas the 3D test on 48 cores required 470.27 core-hours, achieving a speed up factor of ∼195. We see very strong agreement in results between the 2D and 3D results, as seen in Fig. 7 and 8. This agreement might seem trivial since the 2D and 3D grids are closely matched, but the polar transformation significantly alters the flux divergence and metric derivative source terms considered separately. Also, factoring is essential for avoiding strong axis artifacts.

IV. CONCLUSION
We have implemented an axisymmetric evolution of the general relativistic hydrodynamics equations through modification of the local coordinate transformations of a multipatch scheme. Without the appropriate factoring of singular terms from spatial derivatives near the symmetry axis, we find that unphysical errors grow in evolved quantities. Testing of this method, with factoring of singular terms applied, produces results that compare favorably to full 3D simulations at a fraction of the required computational time. Since only minimal modification of the implementation of the evolution equations in 3D was required, this method provides a path for a quick application of axisymmetric evolution to codes that make use of computational domains with local coordinate transformations.
We plan to move forward using this method in order to study the effects of a wide variety of physical parameters on binary post-merger environments that require evolution on secular timescales that we have been unable to explore in the past. Additionally, although our method currently evolves Einstein's equations in 3D using spectral methods, we would also like to extend axisymmetry to the evolution of those equations as well.