We introduce a massively parallelized test-particle based kinetic Monte Carlo code that is capable of modeling the phase space evolution of an arbitrarily sized system that is free to move in and out of the continuum limit. Our code combines advantages of the DSMC and the Point of Closest Approach techniques for solving the collision integral. With that, it achieves high spatial accuracy in simulations of large particle systems while maintaining computational feasibility. Using particle mean free paths which are small with respect to the characteristic length scale of the simulated system, we reproduce hydrodynamic behavior. To demonstrate that our code can retrieve continuum solutions, we perform a test-suite of classic hydrodynamic shock problems consisting of the Sod, the Noh, and the Sedov tests. We find that the results of our simulations which apply millions of test-particles match the analytic solutions well. In addition, we take advantage of the ability of kinetic codes to describe matter out of the continuum regime when applying large particle mean free paths. With that, we study and compare the evolution of shock waves in the hydrodynamic limit and in a regime which is not reachable by hydrodynamic codes.