We investigate the effects of (2+1)d event-by-event fluctuating hydrodynamic backgrounds on the nuclear modification factor and momentum anisotropies of heavy flavor mesons. Using the state-of-the-art D and B mesons modular simulation code (the so-called DAB-mod), updated recently with heavy-light quark coalescence, we perform a systematic comparison of different transport equations, including two energy loss models and a relativistic Langevin model with two drag parametrizations. We present the resulting D$^0$ meson $R_{AA}$, $v_2$ and $v_3$, using the multiparticle cumulant method, in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV and compare them to the latest experimental data. We investigate the $v_2\{4\}/v_2\{2\}$ ratio as a function of centrality for different initial conditions (MCKLN vs. Trento) and different system geometries and sizes (coming from Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV, spherical and prolate Xe-Xe collisions at $\sqrt{s_{\rm NN}}=5.44$ TeV).