Using complex Langevin dynamics we examine the phase structure of complex unitary matrix models and compare the numerical results with analytic results found at large $N$. The actions we consider are manifestly complex, and thus the dominant contribution to the path integral comes from the space of complexified gauge field configuration. For this reason, the eigenvalues of unitary matrix lie off the unit circle and venture out in the complex plane. One example of a complex unitary matrix model, with Polyakov line as the unitary matrix, is an effective description of a QCD at finite density and temperature with $N$ number of colors and $N_f$ number of quark flavors defined on the manifold $S^1 \times S^3$. A distinct feature of this model, the occurrence of a series of Gross-Witten-Wadia transitions, as a function of the quark chemical potential, is reproduced using complex Langevin simulations. We simulate several other observables including Polyakov lines and quark number density, for large $N$ and $N_f$ and found excellent match with the analytic results.