  We discuss the ability of the smoothed particle hydrodynamics (SPH) method
combined with a grid-based solver for the Poisson equation to model mass
accretion onto protostars in gravitationally unstable protostellar discs. We
scrutinize important features of coupling the SPH with grid-based solvers and
numerical issues associated with (1) large number of SPH neighbours and (2)
relation between gravitational softening and hydrodynamic smoothing length.
  We report results of our simulations of razor-thin disc prone to
fragmentation and demonstrate that the algorithm being simple and homogeneous
captures the target physical processes - disc gravitational fragmentation and
accretion of gas onto the protostar caused by inward migration of dense clumps.
  In particular, we obtain two types of accretion bursts: a short-duration one
caused by a quick inward migration of the clump, previously reported in the
literature, and the prolonged one caused by the clump lingering at radial
distances on the order of 15-25~au. The latter is culminated with a sharp
accretion surge caused by the clump ultimately falling on the protostar.
