Heavy flavor probes provide important information about the in-medium properties of the quark-gluon plasma produced in heavy-ion collisions. In this work, we investigate the effects of 2D+1 event-by-event fluctuating hydrodynamic backgrounds on the nuclear suppression factor and momentum anisotropies of heavy flavor mesons and non-photonic electrons. Using the state-of-the-art D and B mesons modular simulation code (called "DAB-MOD"), we perform a systematic comparison of different transport equations in the same background, including a few energy loss models --- with and without energy loss fluctuations --- and a relativistic Langevin model with different drag parametrizations. We present the resulting D and B mesons $R_{AA}$, $v_2$, $v_3$, and $v_4$ as well as multi-particle cumulants, in AuAu collisions at $\sqrt{s_\text{NN}}=200\, \text{GeV}$ and PbPb collisions at $\sqrt{s_\text{NN}}=2.76 \, \text{TeV}$ and $\sqrt{s_\text{NN}}=5.02 \, \text{TeV}$, and compare them to the available experimental data. The $v_2\{4\}/v_2\{2\}$ ratio, which is known to be a powerful probe of the initial conditions and flow fluctuations in the soft sector, is also studied in the context of heavy flavor. We also investigate the correlations between the transverse anisotropies of heavy mesons and all charged particles to better understand how heavy quarks couple to the hydrodynamically expanding quark-gluon plasma. We study the influence that different initial conditions and the implementation of heavy-light quark coalescence has on our results.