  We propose a new method for simulating lattice gauge theories in the presence
of fermions. The method combines flow-based generative models for local gauge
field updates and hierarchical updates of the factorized fermion determinant.
The flow-based generative models are restricted to proposing updates to
gauge-fields within subdomains, thus keeping training times moderate while
increasing the global volume. We apply our method performs to the 2-dimensional
(2D) Schwinger model with $N_f=2$ Wilson Dirac fermions and show that no
critical slowing down is observed in the sampling of topological sectors up to
$\beta=8.45$. Furthermore, we show that fluctuations can be suppressed
exponentially with the distance between active subdomains, allowing us to
achieve acceptance rates of up to $99\%$ for the outer-most accept/reject step
on lattices volumes of up to $V=128\times128$.
