Yes they do. Specifically, the moving puncture method in numerical relativity generally solves for regions inside the event horizon. I asserted (in my answer to the motivating question) that they don't have to use this method. But most currently do.
In "Improved Moving Puncture Gauge Conditions for Compact Binary Evolutions" by Zachariah B. Etienne, John G. Baker, Vasileios Paschalidis, Bernard J. Kelly, and Stuart L. Shapiro in Physical Review D. Rev. D 90, 064032 (2014) DOI:10.1103/PhysRevD.90.064032 arXiv:1404.6523 the authors claim (in 2014) that most simulations with compact objects still use the moving puncture method.
Shapiro has a textbook (Numerical Relativity 2010) that describes the method in a readable and accessible fashion. First on page 432 he describes the puncture method as having a conformal factor for a Schwarzschild example as being factored like $\Psi=\left(1+\frac{M}{2r}\right)f$ and they only solve for $f$ and only let $f$ change in time which requires a singularity at the coordinate $r=0$ regardless of how things actually evolve in time. That older method is then contrasted with the more modern "moving puncture" method.
Then, describing the moving puncture method (the one Shapiro describes in 2014 as still being used in most code simulating compact sources) the authors on page 433 state
the metric is evolved in its entirety. The puncture is allowed to move freely in accord with the gauge conditions, except that care is taken that the singularity never hits an actual grid point. Usually this can be accomplished quite easily. In situations that feature equatorial symmetry, for example, the puncture always remains on the (orbital) plane of symmetry. Using a cell centered finite difference scheme (see e.g. Figure 6.2) no grid points are located on this plane, so that the punctures can never encounter a grid point.
For instance, if you put cubes to cover the region $z\geq 0$ and also put cubes for the region $z\leq 0$ and then always have grid points be the centers of cubes, then no grid points will ever be on the xy plane (that plane is always a face of a cube, and the grid points are never on a face). So the grid points avoid the singularities.
The fact that you have to avoid the grid point being a singularity, means you are inside the event horizon. I argued in my answer to the linked question, that this is not required. A reason to not use the method I described is that then your code needs to be able to evolve matter and also solve Einstein's Equation with source terms (stress energy). If you model eternal black holes, your code only needs to model the vacuum Einstein Equation, $G_{\mu\nu}=0$ and there is no matter to evolve.
But that means since the events before the horizon formed still affect you, your simulation of an astrophysical black hole are technically inaccurate (though maybe in a way swamped by your numerical errors) even before you get to the horizon.
So yes, people put grid points inside the horizon, and yes this is even common. And yes, with appropriate care you can get good results and get them in a reliable way that can be accurate enough to make testable predictions. It isn't necessary to do so.
In fairness, I haven't actually explained the details showing that anyone goes inside the horizon. For instance, if you used a Einstein-Rosen bridge solution for Schwarzschild you can see an initial slice of constant Schwarzschild time doesn't actually cross inside the event horizon, it just touches it, effectively the radius outside the horizon connects to the throat of a (non traversable) wormhole where the white hole and black hole horizons touch and connects the two universes. So you might argue that if only you could stay on these slices you'd never go inside the horizons. But people do often go inside, for technical reasons. It's the choice of slicing spacetime that determines whether you end up going inside, and it's how you select your grid points within a slice that determines whether your singularity hits a grid point. But the point is they try to avoid the singularity hitting the grid point. They don't worry about the event horizon. In fact, if they kept the grids outside the event horizon, they would not have to worry about the singularity, since it is inside.
It is clear that if someone is making an effort to avoid a grid point having a singularity, then they haven't selected a method that systematically and automatically avoids going inside the event horizon. And there is no real reason not to go inside, especially if your data is already technically a bit off before you even get to the horizon. Once you choose to ignore that, why not go in, especially if it helps you introduce less errors in your approximations. You are approximating anyway if you ignore the matter than formed the two black holes.
And in my original answer to the linked question I point out that since the original infalling matter is still affecting you (in an extremely time dilated way), it's always the events prior to the event horizon formation that matter.
