Convergence Rate of NDSolve by Increasing the Spatial Grids

For a parabolic PDE as yours, the a priori error estimates are typically of the form $\mathrm{error} \sim (h^k + \tau^{k/2})$ for some $k \geq 0$ that depends on both the method and the norm in which you measure the error. Here $\tau >0$ is the time step size (i.e. total time divided by number of spatial grids) and $h>0$ is the maximal cell size in the spatial grid. The error goes to $0$ only if both $h \to 0$ and $\tau \to 0$, but you let go only $h \to 0$ while you keep the choice of $\tau$ to NDSolve by using MaxSteps -> Automatic. It it might not converge to $0$.

Edit:

After a some redious reverse engineering, I do not understand what the problem is. I get a nicely decaying sequence of relative $L^2$-errors from the following:

cA[x_, t_] := Erf[x/(2 Sqrt[t])]
pGoal = 8;
xMin = tMin = 0;
xMax = 10;
tMax = 1;
order = 4;

Table[
 xgrid = Subdivide[N@xMin, xMax, nx];
 sol = NDSolveValue[{
    D[u[x, t], t] == D[u[x, t], x, x],
    u[xMax, t] == 1,
    u[x, tMin] == If[x == xMin, 0, 1],
    u[xMin, t] == 0
    },
   u,
   {x, xMin, xMax}, {t, tMin, tMax},
   MaxSteps -> Automatic,
   InterpolationOrder -> Automatic,
   PrecisionGoal -> pGoal,
   Method -> {
     "MethodOfLines",
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "DifferenceOrder" -> order, "Coordinates" -> {xgrid}},
     "DifferentiateBoundaryConditions" -> Automatic}
   ];
 Divide[
  Sqrt[NIntegrate[Abs[sol[x, t] - cA[x, t]]^2, {x, xMin, xMax}, {t, tMin, tMax}]],
  Sqrt[NIntegrate[Abs[cA[x, t]]^2, {x, xMin, xMax}, {t, tMin, tMax}]]
  ],
 {nx, {100, 200, 400, 800}}]

{0.000645067, 0.000229846, 0.0000890372, 0.0000484648}

I wish to add two things, separating the error of the discrete solution computed by NDSolve from the interpolating error between the interpolation grid, and speeding up the computation of the norm. The interpolation error is noticeable, but not significant in the $L^2$ norm, although it is probably the explanation of the OP's original difficulty with the point-wise-relative 1-norm.

One of the problems with speed is cA, which is undefined at t == 0. A different definition, especially with Compile speeds up computation.

ClearAll[cA];
cA = Compile[{{xt, _Real, 1}}, (* call: cA[{x, t}] *)
   Module[{x = First[xt], t = Last[xt]},
    If[x == 0,
     0.,
     If[t == 0,
      1.,
      Erf[x/(2 Sqrt[t])]
      ]]],
   RuntimeAttributes -> {Listable}, Parallelization -> True];

Some parameters. I memoized the solutions so I could play with them without recomputing them. It's unneeded but some of the rest of the code assumes calling sol[nx] won't be slow.

pGoal = 8;
xMin = tMin = 0;
xMax = 10;
tMax = 1;
order = 4;

nxList = {25, 100, 400, 1600, 6400, 25600}; (* discretization sequence *)

ClearAll[sol];
mem : sol[nx_] := With[{xgrid = Subdivide[N@xMin, xMax, nx]},
   mem = NDSolveValue[
     {D[u[x, t], t] == D[u[x, t], x, x],
      u[xMax, t] == 1, u[x, tMin] == If[x == xMin, 0, 1], 
      u[xMin, t] == 0},
     u, {x, xMin, xMax}, {t, tMin, tMax},
     MaxSteps -> Automatic, InterpolationOrder -> Automatic, 
     PrecisionGoal -> pGoal,
     Method -> {"MethodOfLines", 
       "SpatialDiscretization" -> {"TensorProductGrid", 
         "DifferenceOrder" -> order, "Coordinates" -> {xgrid}}, 
       "DifferentiateBoundaryConditions" -> Automatic}]
   ];

Do[sol[nx], {nx, nxList}] (* optional: pre-compute & memoize solutions *)

We compute the integral for the $L^2$ norm from the solution grid used by NDSolve. The values we need, except for "Grid", are stored in the solution and can be obtained from the InterpolatingFunction; the "Grid" is computed efficiently from the "Coordinates". The value of the "Coordinates" has the form xx = {{x0, x1,..., xj}, {t0, t1,..., tk}}, that is, a list of the x-grid and t-grid.

traprule[yy_, xx_] := 
  Fold[#2.MovingAverage[#, 2] &, yy, Differences /@ xx];

Table[With[{
   xt = sol[nx]@"Coordinates",
   exact = cA@ sol[nx]@"Grid",        (* exact values on grid *)
   approx = sol[nx]@"ValuesOnGrid"},  (* computed solution on grid *)
  Divide @@ {
    traprule[(approx - exact)^2, xt] // Sqrt,
    traprule[exact^2, xt] // Sqrt
    }
  ],
 {nx, {25, 100, 400, 1600, 6400, 25600}}]
ListLogPlot[%, Joined -> True]

(*  {0.00202437, 0.000244795, 0.0000493161, 0.0000394941, 0.000039159, 0.0000393847}  *)