  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).
