  A numerical method for simulation of bubble dynamics in three-dimensional
potential flows is presented. The approach is based on the boundary element
method for the Laplace equation accelerated via the fast multipole method
implemented on a heterogeneous CPU/GPU architecture. For mesh stabilization, a
new smoothing technique using a surface filter is presented. This technique
relies on spherical harmonics expansion of surface functions for bubbles
topologically equivalent to a sphere (or Fourier series for toroidal bubbles).
The method is validated by comparisons with solutions available in the
literature and convergence studies for bubbles in acoustic fields. The accuracy
and performance of the algorithm are discussed. It is demonstrated that the
approach enables simulation of dynamics of bubble clusters with thousands of
bubbles and millions of boundary elements on contemporary personal
workstations. The algorithm is scalable and can be extended to larger systems.
