Numerical relativity simulations of binary black hole inspirals and mergers played a
crucial role in the calculations of the gravitational wave signals
that were just observed by the advanced LIGO and Virgo detectors.
For supermassive binary black holes, these signals might also be accompanied by
observable electromagnetic signals from surrounding accreting gas.
I will briefly review the history of simulation efforts to model these
systems in their astrophysical environments. I will also present some exciting
new results in the context of magnetohydrodynamical simulation, indicating
that supermassive binary black hole sources might indeed be also detectable
in the electromagnetic spectrum in the not too distant future.