The basic FPE is
here, with Dirac initial condition for the applications I mentioned.
If a boundary is reachable and the associated parameter is mean-reverting, many finance models will use a zero-flux boundary condition. For example, see my Vol. II Sec 7.9 for discussion with the Heston model -- although that discussion is applied to getting an analytic soln, not fdm.
In general, boundary issues are model-specific, so a comprehensive discussion would be hard to find.
Another example. For the lognormal SABR model, the parameter/coordinate [$]\sigma = 0[$] is not reachable. In that model (pretty sure), numerical absorption at all spatial boundaries will probably work well.
You're right to worry about the Dirac i.c. Personally, if doing the FPE for [$]p(t,\vec{x}_t | \vec{x}_0)[$], I will generally try to place the starting hotspot of the particle [$]\vec{x}_0[$] exactly on a node and then take the i.c. to be a constant there and zero at all other nodes. Then, the constant is an appropriate function of the (local) lattice spacing to imitate the n-dimensional Dirac delta on the fdm lattice.
The challenge is to get an efficient, apparently convergent, norm-preserving, positivity-preserving solution -- not so easy for many standard 2D models -- IMO.