Predictive Variance of a Gaussian Process Suppose $f$ is a function of some variable say $x$ ($x$ could be multi-dim). Then the GP assumption is written as follows
$$f∼GP(m,k)$$
where $m$ is the mean function and $k$ is the covariance function. It is by definition
$$f(x_1),…,f(x_n)\sim N(m(X),k(X,X))$$
where $X:=[x_1,...,x_n]$.
And now, I want to say something about the predictive variance of new points as $n\rightarrow\infty$. 
Under the Gaussian process assumption, I can say that for new points $X^*:=[x^*_1,...,x^*_n]$ we have that
$$f(x_1^*),f(x_2^*),...,f(x_m^*)|f(x_1),f(x_2),...,f(x_n)\sim N(m^*(X^*),k^*(X^*,X^*))$$
And so what I am curious to know, is what happens to my predictive variance as the my sample size $n\rightarrow\infty$?  Does the predictive variance go to 0? Also, whatever the solution is, it would be very beneficial to also see a source for it.
 A: It's going to depend on your covariance kernel $k(s,t)$. Imagine the trivial case where $k(s,t)=\delta(s-t)\sigma^2$, or white noise. And suppose I sample from, WLOG, $[0,1]$. Then no matter how fine my sampling grid, the variance of the predicted value for some $t^*$ will be $\sigma^2$. My prior, noisy observations are uninformative. Taking more of them does not help my case.
Now suppose that $k(s,t)$ is differentiable In that case, realizations of the process will be differentiable also. Increasing the grid so that values get closer to $t^*$ will effectively put a bound on where $x(t^*)$ can reasonably be, thus shrinking the variance. Where a Taylor series expansion would allow you to get close to your predicted outcome at $t^*$, then a probability argument will allow you to hone in on a prediction with probability 1.
Later edit. Let me clarify and respond to some comments. $k(s,t)$ does not tend to 0 as $s$ tends to $t$. It tends towards the marginal variance of the process at $t$. If it's a stationary process, standardized to have unit variance, then $k(t,t)=1$.
Remember that as long as you have a finite sample, the distribution at the sample points will follow a multivariate normal distribution with variance matrix $\Sigma_n$, say. If your sequence of observations tends towards $t^*$, or becomes dense on $[0,1]$, then $\Sigma_n$ tends towards singularity. This is where differentiability of the kernel comes in, because basically, the columns of $\Sigma_n$ are heading towards collinearity as the sampling points get closer together. This is a consequence of Taylor's theorem. You can readily convince yourself with a few simulations that the matrix becomes ill conditioned for fairly small values of $n$, when you have the very smooth gaussian kernel. So this is real.
Basically, what you are doing is interpolating unseen values of the function in light of what you observe.The variance of the interpolant at $t^*$ is $k(t^*,t^*)-K(t^*,T)'\Sigma_n^{-1}K(t^*, T)$, using equation 2.19 in Rasmussen and Williams. $K(t^*,T)$ is the vector you get by evaluating the covariance between $t^*$ and the $n$ observations you already have. However, this formula does not work in the limiting case, since the $\Sigma$ matrix is tending towards singularity, as I already mentioned.
But that doesn't mean that you can't estimate $x(t^*)$. It will be the limit of the $x(t)$ as $t$ tends towards $t^*$. And this will be certain. I can't remember how you prove this, although it is intuitively obvious. If I find a reference, I will edit my answer. I seem to recall that the proof relies on Taylor series.
However, if your sample values are not dense -- say because your support is unbounded or they are not dense in the neighbourhood of the value you want to predict, the predicted variance will not tend towards 0. You will lack certainty about the outcome at $t^*$. The distribution of the your prediction is simply the conditional distribution at $t^*$ given the outcomes you have already observed. 
Imagine that you draw a curve on a board, and you then pound a bunch of nails along the curve. Then pass an elastic band along the nails. The closer the nails get, the less difference there will be between the band and the curve you drew. 
A: tldr; For fixed domain asymptotics and common covariance function the variance drops to zero, while it is not the case for increasing domain asymptotics.
Correct answer depends on the design i.e. on the points you have $\{x_1, \ldots, x_n\}$. 
There are two common type of designs to consider for $n \rightarrow \infty$: fixed domain asympotitc and increasing domain asymptotics. See the comparison of them below (idea of the figure is from the presentation).  
Let us start with the fixed domain asymptotics: all points belong to a fixed bounded domain.
In this case if the design is dense as $n \rightarrow \infty$, then it is easy to see that variance goes to zero if there is no noise and the covariance function is continuous near the origin. Really, we can upper bound the variance by using only the closest point $x'$ to predict the values near $x$. According to the formula for variance in this case we get the variance $k(x, x) - k(x, x') k(x', x')^{-1} k(x', x) \rightarrow 0$ as $|x - x'| \rightarrow 0 $ for a dense sample.
Now let us consider increasing domain asymptotics.
In this case the domain of interest increases as the sample size $n$ increases e.g. $\{x_1 = h \cdot 1, \ldots, x_n = h \cdot n\}$, here $h$ is a grid step size, $h > 0$. For variance in this case we can get positive value even for infinite sample like $\{\ldots, -2h, -h, 0, h, 2h\}$. In this case for the squared exponential covariance function mean variance for points from $[0, 1]$ is of order $h \exp \left(-\frac{1}{h^2} \right)$. So, even for the infinite sample for a not dense enough sample e.g. $h = 1$ the variance doesn't reduce to zero. For a more precise statement see e.g. Corollary 2 in my recent paper from AISTATS
A: Trying to read in between lines, it seems you are interested in a noiseless Gaussian process.  I do not see a valid question here. You have already assumed your covariance function is $k$.  Nothing happens to it as $n \rightarrow \infty$.
P.S. You have mistakes in your last equation.  The prediction is not about $x_*$, it is about $f(x_*)$; $x_*$ is known.
