Modeling the complex physics in numerical simulations requires the solution of a set of partial or ordinary differential equations. My research interests involve the numerical methods and tools to solve those equations. Our developement and work in this field focus on Godunov type codes which include approximate Riemann solver, like the HLLD Riemann solver. In Flock et al. 2010 we adapted the Godunov code PLUTO for 3D global MHD simulations. The reconstruction of the EMF, see figure on the left, becomes very important especially for the hybrid MHD schemes which solve the magnetic field at the cell interfaces, e.g. like the constrained transport (CT) method. As the MRI instability is a subsonic instability it is of great importance to have a low dissipative scheme to accurately resolve the instability. In this work we find a combination of spacial reconstruction, slope limiter and Riemann solver which is accurate and robust at the same time.

In global disk simulations the dominating flow velocity is the Keplerian velocity. In MRI turbulent thin disks the turbulent velocity is of the order of 1% of the Keplerian flow but the strong azimuthal flow will determine the timestep. With the FARGO MHD scheme, the equations are rewritten in such a way that the mean rotation velocity is substracted. Only the important turbulent velocities appear in the equation. With this scheme, published for PLUTO by Mignone et al. 2012, the computational costs are strongly reduced. A useful side-effect of this method is the reduction of intrinsic numerical dissipation and so the resolving of smaller turbulent scales (right plot).

The detailed dynamical evolution of grains embedded in the hydrodynamical simulations are crucial to understand the gas and dust evolution including their feedback. In our work Mignone, Flock et al. 2019 we developed a Lagrangian particle solver to include particles with different stopping times in our simulations. With the new exponential midpoint method we substantially advance the performance of the solver. Simulations of the streaming instability in a global disk are shown on the left.

In accretion disk simulation only with the solution of radiative transfer one is able to investigate correct temperatures in both optical thin and thick regimes. The solution of radiative transfer requires implicit methods as the gas sound speed is much smaller than the speed of light. The inversion of large matrices becomes difficult especially in spherical geometry where the matrix becomes non-symmetric. In our work Flock et al. 2013 we developed a preconditioned BiCGSTAB solver to solve the radiative transfer equations. With this tool we are able to reach high parallel and serial performances which are usually bellow the computational costs of the MHD equations alone.