Numerical simulations of the hydrodynamics and mass transfer involved in the rapid growth of large KDP crystals (linear dimensions up to 60 cm and masses greater than 300 kg) have been performed. The simulations are fully three-dimensional and time-dependent, and the computational geometry corresponds very closely to the 1000-1 crystallizers currently in use at Lawrence Livermore National Laboratory (LLNL) for the growth of KDP crystals for the National Ignition Facility (NIF). The Reynolds number of the flow is of order 10(5), and the Schmidt number is of order 10(3). The flow is turbulent and dominated by flow separation from the rectangular corners of the rotating crystal and a secondary flow driven by Ekman boundary layers on the support platform. The full three-dimensional structure of the flow field is described, but the major emphasis is placed on understanding the surface shear stress distribution on the crystal, since it is critically involved in the process of inclusion formation. The temporal and spatial evolution of the surface shear stress and its effect on the surface morphological stability is compared for several different crystal sizes and rotation conditions.