Abstract The temperature gradients and distribution evolution within the crystal domain in Czochralski crystal growth process have important role in produced crystal's quality. Precise and tight regulation of temperature distribution and gradients is the most promising approach to ensure the crystal quality. In this work, the coupled crystal pulling dynamics and heat transfer models in Czochralski crystal growth with a time-varying boundary is considered. The moving boundary finite element method is used to obtain the optimal reference temperature trajectory associated with the reference crystal shape taking into account constraints on input and the temperature gradients. The obtained reference trajectory is used to implement a model predictive control strategy to track the reference temperature despite uncertainties in crystal domain geometry evolution and disturbances. Finally, the proposed method is implemented on a high fidelity finite element process model with non-planar interface and results are presented to validate the success of the proposed methodology.