Consider the stochastic differential equation
$$ \textrm{(Ito)}\quad dX_t=f(X_t)dt + g(X_t)dw $$
whose probability distribution is described by the partial differential equation
$$ p_t = -(fp)_x+\left(\frac{g^2}{2}p\right)_{xx} $$
Are there any numerical solution techniques for the SDE that are guaranteed to preserve stationary distributions for finite time step lengths?
For instance, if $f = g g'$, then $ p(x, t) = 1 $ is a stationary solution. However, if I use the regular Euler-Maruyama method on a nonlinear $g$, I get artificial clustering around regions where $g''$ is large.