The high demand of effective heat removal from electronic devices has drawn significant interests in exploring ultrahigh thermal conductivity materials and a better understanding of thermal transport across interfaces. We developed a deviational Monte Carlo method to study the phonon transport with a full phonon scattering matrix in time, real, and reciprocal spaces. Our method uses inputs from first-principles calculations and explicitly calculates the spatial variation of phonon distribution function, thus can accurately simulate time-space dependent heat transport in various materials. Graphitic materials have ultrahigh thermal conductivities. The phonon transport in these materials is in the hydrodynamic regime, a new regime with unique thermal transport characteristics that are not possible in better known ballistic and diffusive regimes. The transport phenomena are fluid-like as can be seen in phonon Poiseuille flow, phonon Knudsen minimum, and second sound. We studied the characteristics of phonon Poiseuille flow in suspended graphene by introducing the concept of phonon hydrodynamic viscosity and proposed a decomposition framework to quantify the contribution from each transport regime. Also, we quantitatively predicted the transient propagation of second sound in bulk graphite and observed lattice cooling effect near the adiabatic boundary by pulse heating. Our studies provide fundamental insights on heat transport in ultrahigh thermal conductivity materials and phonon hydrodynamics in graphitic materials. The interfacial transport phenomena have drawn significant interest but mostly been studied in the Landauer framework which neglects internal phonon scattering and non-equilibrium near the interface. This may explain the large discrepancy commonly observed between experimental data and theoretical predictions. The strong non-equilibrium is a result of complex interplay between the interface scattering and internal phonon scattering. This non-equilibrium distribution decays with distance from the interface and recovers to the bulk phonon distribution at a distance which we define as an effective interfacial region. We find that the internal phonon scattering within the interfacial region provides an important contribution to overall interfacial resistance. Our study provides insights into large discrepancies between experimentally measured interfacial resistances and those calculated from the Landauer formula, which are often found in the literature, thus providing a useful way to interpret experimental data.