The effect of non-axisymmetric magnetic perturbations and E × B rotation on neoclassical toroidal viscosity (NTV) is investigated using a drift-kinetic δf Monte-Carlo simulation code, FORTEC-3D, and the simulation is benchmarked with an analytic formula which uses bounce-average approximation. Although the δf code agrees with the analytic formula if the E × B velocity is low or the radial position is away from the resonant rational flux surface, a clear difference appears in the radial profile of NTV when the E × B velocity becomes large. A double-peak profile of NTV appears around the resonant rational flux surface only in the δf simulation. The double peak is created as a result of the resonance of E × B drift and passing particle motion. The benchmark result suggests that the precise drift-kinetic simulation, which treats both trapped and passing particle contributions to neoclassical viscosity, is essential for quantitative evaluation of the rotation damping rate by NTV when the E × B rotation is not slow.