Simulation of the localization and development of plastic shear bands in fluid-saturated rocks is considered using a nonlinear poroelastoplastic model generalizing Biot’s model for a two-phase fluid-saturated porous medium under small and finite strains. A Drucker-Prager yield criterion and a non-associated plastic flow rule are applied to describe an accumulation and localization of plastic strains in a rock. Additionally, a nonlinear dependence of the model parameters (elastic moduli, Biot’s modulus, permeability, etc.) on porosity is considered as well as a dynamic variation of porosity due to the volumetric deformation of the pore space. An isoparametric spectral element method is used to discretize a geometric model and PDEs on curvilinear unstructured meshes of high order in space. A distinctive feature of the developed algorithm for numerical solving the system of nonlinear PDEs of poroelastoplasticity is the use of the dynamic relaxation method, which provides a quasi-stationary solution using an explicit time integration scheme and an optimal choice of the damping parameter. The suggested algorithm allows efficient implementation on a massively parallel high-performance computing system using CUDA technology. The spectral element mesh is naturally mapped onto the CUDA Grid representing GPU’s multiprocessors, and accordingly, each spectral element is mapped onto a streaming block, within which element’s internal nodes are processed by the corresponding threads of the block. Numerical results of solving a series of model problems of the development of plastic shear bands nearby a borehole drilled in a porous fluid-saturated rock are presented. The dynamic variations of porosity and permeability because of the accumulation of plastic deformations are analyzed.