  Internal flows inside gravitationally stable astrophysical objects, such as
the Sun, stars and compact stars are compressed and extremely subsonic. Such
low Mach number flows are usually encountered when studying for example dynamo
action in stars, planets, the hydro-thermodynamics of X-ray bursts on neutron
stars and dwarf novae. Treating such flows is numerically complicated and
challenging task. We aim to present a robust numerical tool that enables
modeling the time-evolution or quasi-stationary of stratified low Mach number
flows under astrophysical conditions. It is argued that astrophysical low Mach
number flows cannot be considered as an asymptotic limit of incompressible
flows, but rather as highly compressed flows with extremely stiff pressure
terms. Unlike the pseudo-pressure in incompressible fluids, a Possion-like
treatment for the pressure would smooth unnecessarily the physically induced
acoustic perturbations, thereby violating the conservation character of the
compressible equations. Moreover, classical dimensional splitting techniques,
such as ADI or Line-Gauss-Seidel methods are found to be unsuited for modeling
compressible flows with low Mach numbers. In this paper we present a nonlinear
Newton-type solver that is based on the defect-correction iteration procedure
and in which the Approximate Factorization Method (AFM) is used as a
preconditioner. This solver is found to be sufficiently robust and is capable
of capturing stationary solutions for viscous rotating flows with Mach number
as small as $\mcal{M} \approx 10^{-3},$ i.e., near the incompressibility limit.
