Complex compressible fluid flows arise in many applications. In the case of inviscid flows, we obtain the Euler equations of gas dynamics that belong to the class of hyperbolic conservation laws. Multidimensional systems of hyperbolic conservation laws pose real challenges both in the analysis as well as in their simulations. Indeed, multidimensional Euler equations are known to be ill-posed in the class of weak entropy solutions since infinitely many solutions for given initial data may exist. Consequently, the convergence of numerical methods is a challenging question. We will introduce the concept of very weak, dissipative solutions, and show that in general well-known structure-preserving numerical methods converge only weakly to these solutions. We also show how to turn weak convergence into strong convergence by applying a set-valued compactness approach due to Komlos. In addition, if a strong solution exists the relative energy can be used to compute error estimates. Another challenging problem is the simulation of compressible flows with uncertain data or model parameters. The latter leads to a random system of partial differential equations. We discuss our recent results on random compressible Navier-Stokes equations and their approximation by Monte Carlo combined with a structure-preserving finite volume method. Cloud simulation is a challenging problem arising in meteorology. We consider random weakly compressible Navier-Stokes equations coupled with a cloud system that models the dynamics of warm clouds. The time evolution of data uncertainties is approximated by the stochastic Galerkin method that combines an asymptotic preserving finite volume method with the spectral-type approximation based on the generalized polynomial chaos expansion in the random space. Extensive numerical experiments confirm second-order accurate approximation in both space and time and exponential convergence in the random space.