Grain boundary migration in the presence of thermal grooving effects play a critical role in the stability of the thin polycrystalline films used in numerous technological applications. We present a computational framework for simulating this motion which relies on a physical model due to Mullins [33, 26], according to which the grain boundaries and the external surfaces are governed by mean curvature motion and surface diffusion, respectively, and along the thermal grooves where the grain boundaries and the exterior surfaces couple, balance of mechanical forces, continuity of the chemical potential, and balance of mass flux dictate the boundary conditions. By adopting an equi-spaced parametric description for the evolving surfaces [40, 13], the physical model can be formulated as a coupled systems of ODEs and PDEs. We propose a finite difference algorithm based on staggered grids for solving the resultant system, which we implement on a polycrystalline layer with columnar structure containing three grains, assuming isotropy. Our algorithm, which is second order accurate in space and first order accurate in time, conserves mass, dissipates energy, and satisfies the predictions of Mullins[33], von Neumann-Mullins [32, 51] and Genin, Mullins, Wynblatt [22], in appropriate limits. Effects such as pitting, hole formation, grain annihilation, wetting, and dewetting can be analyzed using our approach [15, 14, 17, 16].