Abstract It is well-known that magnetohydrodynamics (MHD) dominates the dynamic of the liquid metal flows inside the breeding blankets (BB) of future nuclear fusion plants by magnetic confinement. MHD is a multiphysics phenomenon involving both electromagnetism and incompressible fluid mechanics. From the computational point of view, the simulation of MHD flows in fusion relevant conditions entails a significant challenge. Indeed, due to the shape of the induced electrical currents inside the bulk of the fluid, high spatial resolutions are needed to capture the large gradients found in boundary layers and 3D effects. Besides, solving the equations accurately typically requires very small time steps for the transient algorithms. Over the past few decades, some parallel MHD codes have been developed with success to simulate complex flows in increasingly realistic geometries. Among them, the MHD tools of commercial CFD platforms have attracted attention due to their relatively soft learning curve. Most of these codes are based on the so called ϕ-formulation which, by applying the divergence free condition of the current density to the Ohms law, reduces the electromagnetic part of the problem to a single Poisson equation. As a downside, the approach segregates the fluid and electromagnetic problem. In practice, this establishes important limits to the mesh element size, to the mesh quality and to the time-step needed to obtain accurate and stable solutions that maintains charge conservation at a discrete level. In this work, these limits are explored for the commercial platform ANSYS-Fluent using a test geometry under different conditions. As an alternative, a new code based on Finite Element Methods (FEM) is introduced as well. This open-source code, called GridapMHD (https://github.com/gridapapps/GridapMHD.jl), aims at solving the full set of MHD equations using a monolithic approach. GridapMHD is still in early stages of development but it has already shown promising results.