The famous domino shuffling algorithm was invented to generate the domino tilings of the Aztec Diamond. Using the domino height function, we view the domino shuffling procedure as a discrete-time random height process on the plane. The hydrodynamic limit from an arbitrary continuous profile is deduced to be the unique viscosity solution of a Hamilton-Jacobi equation $u_t+H(u_x)=0$, where the determinant of the Hessian of $H$ is negative everywhere. The proof involves interpolation of the discrete process and analysis of the limiting semigroup of the evolution. In order to identify the limit, we use the theories of dimer models as well as Hamilton-Jacobi equations. It seems that our result is the first example in $d>1$ where such a full hydrodynamic limit with a nonconvex Hamiltonian can be obtained for a discrete system. We also define the shuffling height process for more general periodic dimer models, where we expect similar results to hold.