We investigate the properties of dissipative full discretizations for the advection--diffusion equations associated with models of flow and radiative transport inside stars. The fundamental equation of motion is the fully compressible Navier-Stokes equation, which describes momentum conservation. The model is completed by a continuity equation which ensures conservation of mass, and a total energy equation which describes conservation of the latter. The equations incorporate a radiative source term Q_rad which is determined as the stationary limit of the radiative transfer equation. The equations of hydrodynamics are closed by the equation of state which describes the relation between the thermodynamic quantities. For the spatial discretization of the hyperbolic terms, discretizations of ENO (essentially non-oscillatory) type are implemented. These methods use adaptive stencils which are chosen such as to avoid spurious oscillations in the computed solution. The spatial derivatives are calculated for each direction separately. We derive dissipative space discretizations for the parabolic terms and demonstrate that together with specially adapted implicit total-variation-diminishing (TVD) Runge-Kutta time discretizations with adaptive step-size control this yields reliable and efficient integrators for the underlying high-dimensional nonlinear evolution equations. We demonstrate that fully implicit SDIRK Runge-Kutta methods enable large step-sizes as compared to classical explicit integrators and in conjunction with suitable methods for the associated nonlinear algebraic equations have the potential to significantly reduce the computational effort. In the special parameter regime associated with semiconvection, implicit-explicit methods are demonstrated to provide the most efficient solution methods. The methods we recommend based on our analysis provide advantages in efficiency and accuracy as compared to the classical explicit time integrators.