Nice example! But I think the dimensions D, H, W corresponds to the directions z, y, x, respectively. If this is true, then meshx, meshy, meshz = torch.meshgrid((d, d, d)) should be replaced with meshz, meshy, meshx = torch.meshgrid((d, d, d)). I’m not sure if there is another use case though.
From gridsample documentation:

grid[n, d, h, w] specifies the x, y, z pixel locations for interpolating output[n, :, d, h, w]

So indeed meshx should be put first when stacked. But since meshx changes in the dimension W, it should be the last of torch.meshgrid outputs.
Here is an example motived by yours. The element values indicate their coordinates, e.g. 291 is at (d,h,w) = (2,9,1).
N, C, D, H, W = 1, 1, 3, 10, 10
x = torch.arange(N*C*D*H*W).view(N, C, D, H, W).float()

d = torch.linspace(-1, 1, 2)
h = torch.linspace(-1, 1, 3)
w = torch.linspace(-1, 1, 4)
meshz, meshy, meshx = torch.meshgrid((d, h, w))
grid = torch.stack((meshx, meshy, meshz), 3)
grid = grid.unsqueeze(0) # add batch dim

out = F.grid_sample(x, grid, align_corners=True)
print(x, grid, out)