The principal life-limiting mechanism of Hall thrusters is erosion of the discharge channel walls via ion sputtering. Few sputter yield measurements exist for hexagonal boron nitride (hBN), a common wall material, in the range of ion energies of interest to Hall thrusters, and past modeling efforts have had limited success due to the high computational costs associated with sputter modeling. Presented in this work is an updated molecular dynamics model for the sputtering of hexagonal boron nitride by energetic xenon ions. The model uses graphics processing units to accelerate its computations, allowing better statistics to be achieved when calculating the sputter yields. The model is used to simulate hBN sputtering by ions with kinetic energies of 100 eV and 250 eV and incidence angles of 0° and 45°. It is demonstrated that nitrogen is rapidly depleted from the hBN lattice, leaving a boron-enriched layer and causing the total sputter yield to be strongly dependent on ion fluence. Of the four cases tested, only the 100 eV, 45° incidence case appears to be close to a steady state after more than 35,000 ion impacts. Comparisons to experimental measurements by weight loss and quartz crystal microbalance suggest that the 100 eV, 45° incidence case agrees well with existing data. Differential sputter yields are computed and are shown to be well-described by the modified Zhang function. The species composition of the sputtered particles is analyzed, revealing that the most common sputtered species are atomic boron and diatomic nitrogen. Velocity distribution functions for these two species are computed. The VDF of atomic boron normal to the hBN surface is shown to follow a Sigmund-Thompson distribution. The VDF of diatomic nitrogen exhibits a strong bimodal nature and is better described by a superposition of flux-biased Maxwell-Boltzmann distributions.