In this work, the details of thermal transport in graphite is studied. The aim is to un derstand a striking anisotropy observed in the thermal conductivity of twisted stacks of van-der-Waals thin sheets, as well as to explore the limitations of the relaxation time approximation of the Boltzmann transport equation. To this end, classical molecular dynamics simulations were performed, and analyzed using projection onto phonon modes. Atomic forces were modeled using neuroevolution potentials previ ously trained on results from density functional theory, and gpumd was used to run simulations on graphical processing units. Phonon modes were defined from Second order force constants extracted from the neuroevolution potential using hiphive. The simulations were done both in equilibrium and in variations of non-equilibrium. The equilibrium phonons were modeled using Langevin dynamics to extract their lifetimes and finite-temperature renormalized frequencies. Non-equilibrium states were studied by exciting individual out-of-plane modes to extreme amplitudes be fore running the simulations, and then following the energy as it was distributed from the pumped modes to the rest of the system. Contrary to the exponential decay predicted by the relaxation time approximation, energy was observed to os cillate back and forth between modes when initially excited to a sufficient degree. Furthermore, the system was observed to reach a long-lived quasi-equilibrium state, where a set of coupled modes collectively retained an above-equilibrium energy for an extended period of time. The quasi-equilibrium lifetime was found to depend negatively on the size of the simulated system, but the limit at large sizes remains unknown.