The MH3D++ code has been used to simulate the 3D dynamics of pellet injection. MH3D++ [1] is a modification of the PPPL MH3D resistive full MHD code [2], which replaces the original finite difference / spectral discretization, with an unstructured mesh and finite element [6],[7] / spectral discretization. The unstructured mesh is made with triangular and quadrilateral cells. A sample mesh similar to that used in the computations, but with only 1/4 as many mesh points for clarity, is shown in Fig.1.
The MHD equations are discretized with a zero residual Galerkin approach, using piecewise linear basis functions. The differential operators in the equations become sparse matrices involving integrals of the basis functions and their derivatives.
The finite element unstructured mesh discretization has been incorporated into MH3D++ with an object oriented approach. The unstructured mesh objects generate an unstructured mesh and produce the sparse matrices which implement differential operators including gradient, curl, and divergence, as well as various Poisson solvers based on an Incomplete Cholesky Conjugate Gradient solver. Also included are mesh operations needed for line integration and contour plots.
An important feature of this approach is that most of the MH3D code is retained. This allows direct benchmarking of the two versions against each other. Equilibrium and stability calculations using the two versions have been compared, and agree well.
To represent the effect of
electron free streaming along magnetic field lines, or
large parallel temperature transport,
the MH3D code uses the
``artificial sound" method [8]. In this approach, an artificial
parallel velocity 
is introduced, satisfying
The constants 
are chosen so that the temperature
relaxes rapidly to a near equilibrium (5).
The temperature diffusivity
is small and improves numerical stability.
In the absence of damping and dissipation, the system
(12), (13) conserves the flux tube averaged pressure
as well as an "artificial energy", so that in the presence
of damping the system approaches an equilibrium
(5). The system (12), (13) models
electron free streamimg along field lines, which is the
cause of temperature equilibration.
The following computations were done on a poloidal mesh of slightly more than 2000 grid points, and at least 9 toroidal Fourier harmonics in addition to the n=0 mode. The plasma is assumed bounded by a rigid conducting wall.
An initial equilibrium was prepared, starting with a prescribed
initial, non equilibrium state, and evolving in 2D,
removing kinetic energy, until an equilibrium
is approached. In the following three cases, the equilibrium has
a rotational transform at the magnetic axis 
initial
peak 
and the aspect ratio R/a = 5.
In addition, the pressure is boosted an additional constant
amount with 
This additional constant pressure enhances
the sound speed,
while having no effect on the background magnetic equilibrium.
This makes the total peak
The initial equilibrium
density is constant. The equilibrium flux surfaces are shown
in Fig.2, and the 
profile in Fig.3.
The initial equilibrium was perturbed with a density blob representing
the ionized pellet ablation cloud. In the following, the blob's
peak density is up to 15 times the background density.
The blob is cigar shaped, with a circular cross section in the
poloidal plane, and a dependence on toroidal angle 
proportional
to 
The density perturbation extends
roughly 1/4 of the way around the torus. The model ablation cloud is
considerably less localized than in experiments. This was done
because of constraints on numerical resolution, that should be
relaxed in future simulations. Similarly the density contrast is
less than in experiments. Nonetheless, the effects demonstrated here
are expected to be sufficiently robust to qualitatively account for
experimental observations.
The temperature was initialized adiabatically, as described in
Section II. Before the pellet is injected, the temperature is
given by 
where 
is the magnetic
flux function. After the pellet is inserted, the density is
flux surface averaged, as in (8), and used to calculate
a new temperature, (9). The new pressure is then obtained
from (3). This gives an initial pellet state with a
pressure highly nonuniform on flux surfaces, although the
flux surface averaged pressure is unchanged. This causes a
3D motion towards a new equilibrium state.
Since the profile 
is the same with and without the
pellet, there is a tendency to return to the initial
equilibrium state after the pellet cloud spreads out on magnetic
field lines. However, when reconnection occurs, the profiles
change and the pellet does not spread out to the initial
equilibrium.
Simulations were also done with a simpler initialization
in which the temperature was not modified by the 
adiabatic correction. These simulations were quite similar in
behavior, as far as the maximum displacement of the pellet cloud
is concerned. However, since in this case the initial
profile was changed, the pellet had less tendency
to return to the initial equilibrium.
