Polydisperse bubbly flows appear in many industrial applications and are particularly relevant for nuclear and process engineering. A popular approach for their simulation is the two-fluid model, which requires information about the bubble size. An associated challenge is to incorporate the effects of coalescence and breakup on the size distribution, which can be tracked by a population balance equation. It can be solved by the method of classes, which splits the bubble population into a finite number of size groups. However, this technique is expensive because the source term assembly involves the computation of coalescence and breakup frequencies between all bubble size pairs. This work focuses on improving the performance of the class method implementation within the OpenFOAM Foundation release (www.openfoam.org), an open source computational fluid dynamics software that allows for domain decomposition and parallel execution on multiple central processing units. Here, the parallel computation of the coalescence and breakup frequencies is shifted to graphics processing units using the Nvidia CUDA framework. Calculations are done asynchronously with overlapping data transfers, treating size pairs as independent units of work that are computed in parallel. The coalescence and breakup frequency computation on graphics processing units leads to a significant speedup compared to the pre-existing implementation. The improvement is demonstrated for a co-current two-phase flow in a vertical pipe. Both the source code and the case setup are made publicly available.