Numerical relativity - coordinates, discretization and shift vector I would like to understand how, in the numerical implementation of the relativistic 3+1 decomposition, the spatial coordinates are evolved from slice to slice and what is their relation to the shift vector.
Let me elaborate on the problem. In 3+1 decomposition, the evolution of spacetime $\mathcal{M}$ can be understood as evolving spatial fields on a hypersurface $\Sigma$ that changes its topology and geometry. In other words, fields and coordinates are evolved simultaneously. 
In the 3+1 framework, there exists a concept of a Eulerian observer.The tangent vector field to Eulerian observers' worldlines is the vector field $\boldsymbol{n}$ normal to $\Sigma$.
Then, there is a concept of a coordinate time observer congruence, whose tangent vector field is $$\boldsymbol{\partial_{t}} = \alpha \boldsymbol{n} + \boldsymbol{\beta} \equiv \boldsymbol{m} + \boldsymbol{\beta},$$ where $\alpha$ is the lapse and $\boldsymbol{\beta}$ is the shift vector.  
The worldlines of coordinate time observers are curves along which the spatial coordinate values are constant. If the curve $\gamma$ corresponds to $\boldsymbol{\partial_{t}}$, then solving the
$\dot{\gamma}=\boldsymbol{\partial_{t}}$ yields a family of curves $\gamma(\tau)=(t(\tau),x^{i}_{0})$, where the spatial coordinates remain constant.
If I am correct, the coordinate values of the Eulerian observers spatial position can be recovered (by virtue of $\boldsymbol{\partial_{t}} =\boldsymbol{m} + \boldsymbol{\beta}$) via $x^{i}_{\boldsymbol{m}}(s)= x^{i}_{0} - x^{i}_{\beta}(s)$.
Here is when I have a conceptual problem and I would like to know if my thinking is correct.
Let us assume we want to discretize the hypersurface, i.e. represent our hypersurface as a grid.
I will use a three dimensional box and for each point use indices $p=(i,j,k)$.

*

*Now, what is the relation between these indices and values of some coordinates?
Do we assign to grid indices the constant coordinate values? In other words, is it the case that $x^{i}(p)=(i,j,k) = x^{i}_{0}$ and remain constant through the computations?

*I suspect that what could be done is assigning the initial spatial coordinate values to grid points (i,j,k), setting up the initial data (where the spatial initial coordinates of Eulerian and coordinate time observers coincide) and then, at each step in time ($t\rightarrow t + \delta t)$, updating the spatial coordinate values of the grid through the shift vector as prescribed by $x^{i}_{\boldsymbol{m}}(s)= x^{i}_{0} - x^{i}_{\beta}(s)$. Then we would (of course?) have to assume that ANY other quantities are represented in the code in such a way, that they depend on the changing "background" coordinate grid.

*Otherwise, we can assume that the numerical grid positions (i,j,k) always correspond to points that are occupied by coordinate time observers, therefore they never change, as by definition such observers have constant spatial coordinates. 
What is the role played by the shift vector in the simulation then, except for influencing the evolution of other quantities (spatial metric $\gamma_{ij}$, extrinsic curvature $K_{ij}$)? 
Surely not changing the coordinate values of grid points anymore?
The numerical grid stays constant, but it represents the evolving spacetime intrinsically in a sense, implicitly, through the constant spatial coordinate values of the moving (coordinate time) observers.

* If question 2 has it correct, are the values of spatial coordinates on a grid position (i,j,k) monitored in some way by following the coordinate time observers? 
It is difficult to express my question in a consistent way, so I hope it is understandable. I would be grateful for a detailed explanation of how it is done in practice and if my ideas are consistent.
Here is an illustration:
https://imgur.com/a/rAfH0f2
EDIT: I have decided to cross out points 2 and 4, as they seem not to be the case, based on answers in the post, the provided literature and after some consideration.
 A: The short answer is that this is one of the free choices that the numerical relativist has to make when constructing the simulation.  There is no a-priori correct answer as the lapse and shift are exactly related to your free choice of coordinates in general relativity.
For a bit more detail and context:
In the 3+1 decomposition, the lapse and shift are also fields that evolve.  Unlike the physical fields, you can give them any evolution equation that you'd like, including the trivial equation that they don't evolve, but in principle at least they change with coordinate time.
A lot of the early work in numerical relativity was focused on choosing coordinate conditions that would allow for stable numerical solutions, especially in black hole spacetimes.  Without getting into details of which worked and which didn't - a whole field onto itself - imagine if you were starting at the beginning.  Do you want to take a lapse that goes to 0 near the singularity so that you don't have to deal with an infinity in your numerical grid?  That might help with the singularity but your slices will stretch as the evolution proceeds.  Do you want to keep the coordinates from "falling" into the black hole by using a big shift vector?  If you have two black holes, do you want to use coordinates that keep the black holes fixed in your grid ("co-rotating") or let he black holes move to have a simpler shift?  Whatever you want to do, how do you match it up to the coordinates used in the initial data?
Also keep in mind that there are really two problems in numerical relativity.  The one you've focused on is to evolve the fields given initial data.  The other, is to generate the initial data in the first place, which involves solving elliptic constraint equations on the initial slice.  That will introduce its own issues with coordinates for you to ponder, starting with how you understand how that initial slice fits into the full spacetime.
