Since the death of Prof. Dr. Jan Czochralski nearly 50 years ago, crystals grown by the Czochralski method have increased remarkably in size and perfection, resulting today in the industrial production of silicon crystals about 30 cm in diameter and two meters in length. The Czochralski method is of great technological and economic importance for semiconductors and optical crystals. Over this same time period, there have been equally dramatic improvements in our theoretical understanding of crystal growth morphology. Today we can compute complex crystal growth shapes from robust models that reproduce most of the features and phenomena observed experimentally. We should care about this because it is likely to result in the development of powerful and economical design tools to enable future progress. Crystal growth morphology results from an interplay of crystallographic anisotropy and growth kinetics by means of interfacial processes and long-range transport. The equilibrium shape of a crystal results from minimizing its anisotropic surface free energy under the constraint of constant volume; it is given by the classical Wulff construction but can also be represented by an analytical formula based on the ξ-vector formalism of Hoffman and Cahn. We now have analytic criteria for missing orientations (sharp corners or edges) on the equilibrium shape, both in two (classical) and three (new) dimensions. Crystals that grow under the control of interfacial kinetic processes tend asymptotically toward a "kinetic Wulff shape," the analogue of the Wulff shape, except it is based on the anisotropic interfacial kinetic coefficient. If it were not for long range transport, crystals would presumably nucleate with their equilibrium shape and then evolve toward their "kinetic Wulff shape." Allowing for long range transport leads to morphological instabilities on the scale of the geometric mean of a transport length (typically a diffusivity divided by the growth speed) and a capillary length (of the order of atomic dimensions). Resulting crystal growth shapes can be cellular or dendritic, but can also exhibit corners and facets related to the underlying crystallographic anisotropy. Within the last decade, powerful phase field models, based on a diffuse interface, have been used to treat simultaneously all of the above phenomena. Computed morphologies can exhibit cells, dendrites and facets, and the geometry of isotherms and isoconcentrates can also be determined. Results of such computations are illustrated in both two and three dimensions.