Two interacting systems of partial differential equations governing three-dimensional laminar flow of an incompressible viscous fluid undergoing solidification or melting while under the influence of an externally applied magnetic field have been formulated and integrated numerically. The model includes effects of Joule heating, latent heat, and arbitrary magnitude and orientation of gravity and the magnetic field. It allows for arbitrary temperature-dependent physical properties within the melt and the solid phase. The mushy region is captured by varying viscosity orders of magnitude in the mushy region and by allowing latent heat of phase change to be an arbitrary function of temperature. The uniqueness of this approach is in the fact that both liquid and solid phases are treated as incompressible liquids with the solid phase having an extremely high viscosity. It was found numerically that the magnetic field strength and orientation can significantly influence flow field velocity, and vorticity, amount of accrued solid, and the solid/liquid interface shape.