In this paper, we fully characterize the nonlinear evolution of a 2D crystal growing in a supercooled melt using a boundary integral scheme in which a novel time and space rescaling is implemented. Our simulations reveal that at long times there is nonlinear stabilization even though unstable growth may be significant at early times. This stabilization leads to the existence of universal limiting shapes that depend on the applied far-field flux. A number of limiting shapes exist for each flux, but only one is universal in the sense that a crystal with an arbitrary initial shape will evolve to this universal shape. We construct a phase diagram that reveals the relationship between the applied flux and the achievable symmetries of the limiting shapes. The phase diagram can be used to design a nonlinear protocol that might, be used in a physical experiment to control the nonlinear morphological evolution of a growing crystal. Because our analysis shows that interactions among the perturbation modes are similar in both 2D and 3D, our results apply qualitatively to 3D as our preliminary 3D simulation results suggest.