Unexpected 'switch' of partial derivatives when differentiating InterpolatingFunction

Note an answer; just too long for a comment

@user21 's suggestion (taking t,x,y instead of x,y,t as variables) works.
Here is the code :

Ω = Rectangle[{0, 0}, {1, 2}];
tmin = 0; tmax = 10;
Ω = Rectangle[{0, 0}, {1, 2}];
f = NDSolveValue[{
   Laplacian[T[t, x, y], {x, y}] == D[T[t, x, y], t]
   , T[0, x, y] == 1
   , DirichletCondition[T[t, x, y] == y/2, True]
   }
  , T
  , {x, y} ∈ Ω
  , {t, tmin, tmax}
  (*, Method\[Rule]{"MethodOfLines","TemporalVariable"\[Rule] t,
  "SpatialDiscretization"\[Rule]{"FiniteElement"}}*)]

Plot3D[f[tmax, x, y], {x, 0, 1}, {y, 0, 2}, PlotRange -> {0, 1},
 AxesLabel -> {"X", "Y"}]

Echo[Derivative[0, 1, 0][f][tmax, 0.5, 1], 
  "derivative along X axis : "];

Internally, the FEM-generated interpolating function always stores t first. This means that an input where we have an independent variable ordering of x, y, t is internally reordered to t,x,y. There is a chance that something goes wrong during the derivative computation. To try this, change your ordering from x,y,t to t,x,y and see if this fixes the issue. Should that indeed work, then the issue you reported is a bug.

Update:

OK, the proper workaround for this is to put t first as in t,x,y. However I'd like to show another way to check this:

Compute the solution:

tmin = 0; tmax = 10;
Ω = Rectangle[{0, 0}, {1, 2}];
f = NDSolveValue[{Laplacian[T[x, y, t], {x, y}] == D[T[x, y, t], t], 
    T[x, y, 0] == 1, DirichletCondition[T[x, y, t] == y/2, True]}, 
   T, {x, y} ∈ Ω, {t, tmin, tmax}];

Compute and evaluate the derivative:

fx = Derivative[1, 0, 0][f];
fx[0.5, 1, tmax]
0.4999678532261518`

Inspect the interpolating data structure, as explained here:

fx[[2]]
{5, 12417, 1, {63, 1394, 0}, {4, 3, 3}, {0, 0, 1}, 0, 0, 0, 
 Indeterminate &, {}, {}, False}

Note the 6th position:

fx[[2, 6]]
{0, 0, 1}

This gives the derivative order that the interpolating function represents. When we know that the internal representation is t,x,y then {0,0,1} means fy and not fx. Let's change that:

fx[[2, 6]] = {0, 1, 0}
{0, 1, 0}

And....

fx[0.5, 1, tmax]
1.4206200904431552`*^-14

But again: the proper workaround is to use the t,x,y ordering and not the x,y,t one. I found the issue in the code and fixed it. But I am hesitant to include this for 12.1 since the release is probably close and even though all tests pass fiddling with this so short of the release makes me nervous. Should there be a 12.1.1 (and that is not clear at this point in time) I'll consider adding it there. Else it will have to wait for 12.2. Since this has an easy fix I think it's not the end of the world to not have it in 12.1. However, I apologize for my mistake. Sorry.

Update: Should there be a 12.1.1 then the fix for this issue will be in it.