Rejection sampling doesn’t have to be a problem.
Let f\left(y\right) be the probability that the simulated dataset y is rejected (usually a 0-1 function if you have a clear idea what a “bad” dataset looks like, but could be probabilistic if you’re relying on finicky diagnostics). The important numbers are the probability of rejection for parameter \theta
L\left(\theta\right)=\int f\left(y\right)\pi\left(y|\theta\right)\mathrm{d}y
and the total rate of rejections from the prior
R=\iint f\left(y\right)\pi\left(y|\theta\right)\pi\left(\theta\right)\mathrm{d}y\mathrm{d}\theta=\int L\left(\theta\right)\pi\left(\theta\right)\mathrm{d}\theta
Rejecting the parameter draw when it generates a “bad” dataset effectively distorts the prior
\pi\left(\theta\right)\to\frac{L\left(\theta\right)}{R}\pi\left(\theta\right)
and of course rejections change the generating distribution
\pi\left(y|\theta\right)\to\frac{f\left(y\right)}{L\left(\theta\right)}\pi\left(y|\theta\right)
but crucially these changes cancel out when computing the posterior
\pi\left(\theta|y\right)\to\pi\left(\theta|y\right)
Despite the rejections the fits target the correct posteriors and you can still expect the rank histogram to be uniform.