We consider the numerical solution of a phase field model for polycrystallization in the solidification of binary mixtures in a domain $$ \varOmega \subset \mathbb {R}^2$$ . The model is based on a free energy in terms of three order parameters: the local orientation $$\varTheta $$ of the crystals, the local crystallinity $$\phi $$ , and the concentration c of one of the components of the binary mixture. The equations of motion are given by an initial-boundary value problem for a coupled system of partial differential equations consisting of a regularized second order total variation flow in $$ \varTheta $$ , an $$L^2$$ gradient flow in $$\phi $$ , and a $$W^{1,2}(\varOmega )^*$$ gradient flow in c. Based on an implicit discretization in time by the backward Euler scheme, we suggest a splitting method such that the three semidiscretized equations can be solved separately and prove existence of a solution. As far as the discretization in space is concerned, the fourth order Cahn–Hilliard type equation in c is taken care of by a $$\hbox {C}^0$$ Interior Penalty Discontinuous Galerkin approximation which has the advantage that the same finite element space can be used as well for the spatial discretization of the equations in $$ \varTheta $$ and $$ \phi $$ . The fully discretized equations represent parameter dependent nonlinear algebraic systems with the discrete time as a parameter. They are solved by a predictor corrector continuation strategy featuring an adaptive choice of the time-step. Numerical results illustrate the performance of the suggested numerical method.