We propose an adaptive optimization algorithm for operating district heating networks in a stationary regime. The behavior of hot water flow in the pipe network is modeled using the incompressible Euler equations and a suitably chosen energy equation. By applying different simplifications to these equations, we derive a catalog of models. Our algorithm is based on this catalog and adaptively controls where in the network which model is used. Moreover, the granularity of the applied discretization is controlled in a similar adaptive manner. By doing so, we are able to obtain optimal solutions at low computational costs that satisfy a prescribed tolerance w.r.t. the most accurate modeling level. To adaptively control the switching between different levels and the adaptation of the discretization grids, we derive error measure formulas and a posteriori error measure estimators. Under reasonable assumptions we prove that the adaptive algorithm terminates after finitely many iterations. Our numerical results show that the algorithm is able to produce solutions for problem instances that have not been solvable before.