A numerical method for the phase-resolving simulation of rigid particles and deformable bubbles in viscous, magnetohydrodynamic flows is presented. The approach features solid robustness and high numerical efficiency. The implementation is three-dimensional and fully parallel suiting the needs of modern high-performance computing. In addition to the steps towards magnetohydrodynamics, the thesis covers method development with respect to the immersed boundary method which can be summarized in simple words by ¿From rigid spherical particles to deformable bubbles of complex shape¿. Bubble interaction is addressed as well including bubble collision and coalescence and a new coalescence model is introduced. Applications of the method comprise simulations of the rise of a single bubble and a bubble chain in liquid metal with and without a magnetic field highlighting the major effects of the field on the bubble dynamics and the flow field. The effect of bubble coalescence is quantified for two closely adjacent bubble chains.