  In this work, we report a massively parallel and time domain implementation
of the 3D phase field model that can reach beyond micron scale and consider for
arbitrary electrical and mechanical boundary conditions. The first part of the
paper describes the theory and the numerical implementation of the model. A
mixed-mode approach of finite difference and finite element grid has been used
for calculating the nonlocal electrostatic and elastic interactions
respectively. All the local and non-local interactions are shown to scale
linearly up to thousands of processors. This massive parallelization allows us
to directly compare our results with multiple experiments at the same size
scale. The second part of the paper presents results of ferroelectric switching
in devices based on the multi-ferroic BiFeO$_{3}$. We have particularly
emphasized the importance of charge driven domain growth and the effect of
electrical boundary conditions that explain the temporal evolution of
ferroelectric domains observed in recent experiments. We also predict a
mechanism of controlling domain size in the multi-domain ferroelectric
switching that could be useful for practical applications.
