Fusion reactors designs frequently involve the use of liquid metal flows in the presence of strong magnetic fields. Simulation of the flows involves the solution of continuum equations for fluid flow and magnetic induction, usually done with finite difference methods. In this paper, an alternative method, based on the generalized lattice Boltzmann equation (GLBE), and implemented in the MetaFlow code is discussed. It has a number of desirable features, including fast execution, excellent parallel scalability, and can easily handle complex geometries. The use of the recent GLBE variant greatly enhances stability and accuracy. To simulate magnetohydrodynamic (MHD) flows relevant to fusion applications using GLBE, several new models have been developed, including new boundary condition formulations, preconditioners for faster steady-state convergence, variable electrical conductivity materials, and to resolve thin Hartmann layers. These models are discussed, and validations against MHD benchmarks, including 3-D driven cavity, high Hartmann number and turbulent cases are presented.