Internal flows inside gravitationally stable astrophysical objects, such as the Sun, normal and compact stars, are rotating, highly compressed and extremely subsonic. Such low Mach number flows are usually encountered when studying, for example, the dynamo action in stars and planets or the nuclear burst on neutron stars and white dwarfs. Handling of such flows numerically on time-scales longer than the dynamical one is complicated and challenging. The aim of this paper is to address the numerical problems associated with the modelling of internal quasi-stationary, rotating low Mach number flows in stars and to discuss possible solution scenarios. It is shown that the quasi-symmetric approximate factorization method (AFM) as a pre-conditioner within a non-linear Newton-type defect-correction solution procedure is best suited for modelling quasi-stationary weakly compressible flows with moderate low Mach numbers. This method is robust as it can be applied to model time-dependent compressible flows without further modifications. The AFM-pre-conditioning techniques are shown to be extendable into three dimensions with an arbitrary equation of state. Classical dimensional splitting techniques, however, such as the alternating direction implicit or line-Gauss-Seidel methods are not suited for modelling compressible low Mach number flows. It is also argued that hot and low Mach number astrophysical flows cannot be considered as an asymptotic limit of incompressible flows, but rather as highly compressed flows with extremely stiff pressure terms. We show that, unlike the pseudo-pressure in incompressible fluids, a Poisson-like treatment for the pressure would smooth unnecessarily physically induced acoustic perturbations, thereby violating the conservation of the total energy. Results of several hydrodynamical calculations are presented, which demonstrate the capability of the solver to search for solutions, that correspond to stationary, viscous and rotating flows with a Mach number as small as as well as to fluid flows that are subject to ultra-strong Newtonian and general relativistic gravitational field
