I’ve been using Stan to sample from a highly multimodal posterior for a model with 11935 parameters. With my 3+ years old laptop, it takes 2 hours for 1000 warmup and 3000 post-warmup iterations per chain. With GPU this could be get down to less than 40 mins. With some optimization of the GPU code (which is now optimized for big n, and not for big number of parameters) it could be made even faster. The sampling passes all convergence diagnostics (using 4 chains) and results are consistent in repeated runs. Looking at the marginals I estimate there is less than 2^9 modes. Here’s a figure showing the multimodality
The model I used is logistic regression with regularized horseshoe prior. I can email you data and code if you want to test.