We study in this paper three variants of the high-order Discontinuous Galerkin (DG) method with Runge-Kutta (RK) time integration for the induction equation, analysing their ability to preserve the divergence free constraint of the magnetic field. To quantify divergence errors, we use a norm based on both a surface term, measuring global divergence errors, and a volume term, measuring local divergence errors. This leads us to design a new, arbitrary high-order numerical scheme for the induction equation in multiple space dimensions, based on a modification of the Spectral Difference (SD) method [1] with ADER time integration [2]. It appears as a natural extension of the Constrained Transport (CT) method. We show that it preserves $\nabla\cdot\vec{B}=0$ exactly by construction, both in a local and a global sense. We compare our new method to the 3 RKDG variants and show that the magnetic energy evolution and the solution maps of our new SD-ADER scheme are qualitatively similar to the RKDG variant with divergence cleaning, but without the need for an additional equation and an extra variable to control the divergence errors. [1] Liu Y., Vinokur M., Wang Z.J. (2006) Discontinuous Spectral Difference Method for Conservation Laws on Unstructured Grids. In: Groth C., Zingg D.W. (eds) Computational Fluid Dynamics 2004. Springer, Berlin, Heidelberg [2] Dumbser M., Castro M., Par\'es C., Toro E.F (2009) ADER schemes on unstructured meshes for nonconservative hyperbolic systems: Applications to geophysical flows. In: Computers & Fluids, Volume 38, Issue 9

Recent work has suggested that the net gravitational force acting on a massive and luminous perturber travelling through a gaseous and opaque medium can have same direction as the perturber's motion (an effect sometimes called negative dynamical friction). Analytic results were obtained using a linear analysis and were later confirmed by means of non-linear numerical simulations which did not resolve the flow within the Bondi sphere of the perturber, hence effectively restricted to weakly perturbed regions of the flow (paper~I). Here we present high resolution simulations, using either 3D Cartesian or 2D cylindrical meshes that resolve the flow within the Bondi sphere. We perform a systematic study of the force as a function of the perturber's mass and luminosity, in the subsonic regime. We find that perturbers with mass $M$ smaller than a few $M_c\sim \chi c_s/G$ are subjected to a thermal force with a magnitude in good agreement with linear theory ($\chi$ being the thermal diffusivity of the medium, $c_s$ the adiabatic sound speed and $G$ the gravitational constant), while for larger masses, the thermal forces are only a fraction of the linear estimate that decays as $M^{-1}$. Our analysis confirms the possibility of negative friction (hence a propulsion) on sufficiently luminous, low-mass embryos embedded in protoplanetary discs. Finally, we give an approximate expression of the total force at low Mach number, valid both for sub-critical ($M<M_c$) and super-critical ($M>M_c$) perturbers.

We present a two-dimensional Cartesian code based on high order discontinuous Galerkin methods, implemented to run in parallel over multiple GPUs. A simple planet-disc setup is used to compare the behaviour of our code against the behaviour found using the FARGO3D code with a polar mesh. We make use of the time dependence of the torque exerted by the disc on the planet as a mean to quantify the numerical viscosity of the code. We find that the numerical viscosity of the Keplerian flow can be as low as a few $10^{-8}r^2\Omega$, $r$ and $\Omega$ being respectively the local orbital radius and frequency, for fifth order schemes and resolution of $\sim 10^{-2}r$. Although for a single disc problem a solution of low numerical viscosity can be obtained at lower computational cost with FARGO3D (which is nearly an order of magnitude faster than a fifth order method), discontinuous Galerkin methods appear promising to obtain solutions of low numerical viscosity in more complex situations where the flow cannot be captured on a polar or spherical mesh concentric with the disc.