We have developed a numerical algorithm for simulation of wave propagation in double-porosity media, where the pore space is saturated with a single fluid. Spherical inclusions embedded in a background medium oscillate to yield attenuation by mode conversion from fast P-wave energy to slow P-wave energy (mesoscopic or wave-induced fluid-flow loss). The theory is based on the Biot theory of poroelasticity and the Rayleigh model of bubble oscillations. The differential equation of the Biot-Rayleigh variable is approximated with the Zener mechanical model, which results in a memory-variable viscoelastic equation. These approximations are required to model mesoscopic losses arising from conversion of the fast P-wave energy to slow diffusive modes. The model predicts a relaxation peak in the seismic band, depending on the diameter of the patches, to model the attenuation level observed in rocks. The wavefield is obtained with a grid method based on the Fourier differential operator and a second-order time-integration algorithm. Because the presence of two slow quasistatic modes makes the differential equations stiff, a time-splitting integration algorithm is used to solve the stiff part analytically. The modeling has spectral accuracy in the calculation of the spatial derivatives.