  Hydromagnetic turbulence produced during phase transitions in the early
universe can be a powerful source of stochastic gravitational waves (GWs). GWs
can be modelled by the linearised spatial part of the Einstein equations
sourced by the Reynolds and Maxwell stresses. We have implemented two different
GW solvers into the {\sc Pencil Code} -- a code which uses a third order
timestep and sixth order finite differences. Using direct numerical integration
of the GW equations, we study the appearance of a numerical degradation of the
GW amplitude at the highest wavenumbers, which depends on the length of the
timestep -- even when the Courant--Friedrichs--Lewy condition is ten times
below the stability limit. This degradation leads to a numerical error, which
is found to scale with the third power of the timestep. A similar degradation
is not seen in the magnetic and velocity fields. To mitigate numerical
degradation effects, we alternatively use the exact solution of the GW
equations under the assumption that the source is constant between subsequent
timesteps. This allows us to use a much longer timestep, which cuts the
computational cost by a factor of about ten.
