Consider a stochastic differential equation dY(t)=f(X(t))dX(t), where X(t) is a pure jump Levy process with finite p-variation, 1p. Following the geometric solution construction of Levy-driven stochastic differential equations, we develop a class of epsilon-strong simulation algorithms that allows us to construct a probability space, supporting both Y and a fully simulatable process Y_epsilon, such that Y_epsilon is within epsilon distance from Y under the Skorokhod J1 topology on compact time intervals with probability 1. Moreover, the user can adaptively refine the accuracy levels. This tolerance-enforcement feature allows us to easily combine our algorithm with multilevel Monte Carlo method for efficient estimation of expectations, and adding as a benefit a straightforward analysis of the rate of convergence.