Consider dY(t)=f(X(t))dX(t), where X(t) is a pure jump Levy process with finite p-variation norm, 1<= p < 2, and f is a Lipchitz continuous function. Following the geometric solution construction of Levy-driven stochastic differential equations in (Williams 2001), we develop a class of epsilon-strong simulation algorithms that allows us to construct a probability space, supporting both the geometric solution Y and a fully simulatable process Y_epsilon, such that Y_epsilon is within epsilon distance from Y under the uniform metric on compact time intervals with probability 1. Moreover, the users can adaptively choose epsilon’ < epsilon, so that Y_epsilon’ can be constructed conditional on Y_epsilon. This tolerance-enforcement feature allows us to easily combine our algorithm with Multilevel Monte Carlo for efficient estimation of expectations, and adding as a benefit a straightforward analysis of rates of convergence. This is joint with Jose Blanchet, Fei He and Offer Kella.