Reduce the dimension of InterpolatingFunction

To avoid the (reported) crash you can swap the order of t and the spatial variables like so:

region = Polygon[{{0, 0}, {0, 1}, {1, 0}}];
fksol = NDSolveValue[{Derivative[0, 2, 0][u][t, x, y] + 
      Derivative[0, 0, 2][u][t, x, y] + 
      Derivative[0, 1, 0][u][t, x, y] + 
      Derivative[0, 0, 1][u][t, x, y] == 
     Derivative[1, 0, 0][u][t, x, y], 
    u[0, x, y] == (Erf[x/.1] - Erf[(x - 1)/.1] - 1) (Erf[y/.1] - 
        Erf[(y - 1)/.1] - 
        1) (PDF[NormalDistribution[.2, .1], x]*
         PDF[NormalDistribution[.8, .1], y] // Evaluate), 
    u[t, 0, y] == 0, u[t, x, 0] == 0}, 
   u, {t, 0, 1}, {x, y} \[Element] region];
fksol["ValuesOnGrid"]

Too long for a comment, and possibly an answer....

Question/comment, based on the OP and comments:

The following integrates fksol[0.1, y, t] over the slice of the interpolating function's domain where x == 0.1. Is that what you're after?

Needs["NDSolve`FEM`"]

emesh = fksol["ElementMesh"];

getY[x0_] := 
 With[{emx = ToElementMesh@ DiscretizeRegion@
      RegionIntersection[MeshRegion@emesh, ImplicitRegion[x >= x0, {x, y}]]},
  With[{coords = emx["Coordinates"]},
   With[{bdy = emx["BoundaryElements"] /. LineElement -> List /. 
       idcs : {__Integer} :> coords[[idcs]]},
    Union @@ 
     Cases[bdy, seg : {{x_ /; x == x0, _} ..} :> seg[[All, 2]], Infinity]
    ]]];

tpts = Last@fksol["Coordinates"];

int[x_?NumericQ, y_?NumericQ, t_?NumericQ] := fksol[x, y, t];
With[{x = 0.1},
 NIntegrate[int[x, y, t],
  Evaluate@{y, Sequence @@ getY[x]},
  Evaluate@{t, Sequence @@ tpts}, AccuracyGoal -> 17]
 ]
(*  0.039357  *)