In this paper heavy quark energy loss models are embedded in full event-by-event viscous hydrodynamic simulations to investigate the nuclear suppression factor and azimuthal anisotropy of D$^0$ mesons in PbPb collisions at 5.02 TeV in the $p_T$ range 8-40 GeV. In our model calculations, the $R_text{AA}$ of D$^0$ mesons is consistent with experimental data from the CMS experiment. We present the first calculations of heavy flavor cumulants $v_2{2}$ and $v_3{2}$ (and also discuss $v_2{4}$), which is also consistent with experimental data. Event-shape engineering techniques are used to compute the event-by-event correlation between the soft hadron $v_n$ and the heavy meson $v_n$. We predict a linear correlation between these observables on an event-by-event basis.