I think you want to give more context than that. I assume that phi is your “ground truth loss” (x being the ground truth or some such information) and psi the consistency loss.
Now psi is rather underspecified in the formula. From your description, I looks like that i,j iterates over all pairs of adjacent pixels (of which there would be 8 if you have axis-parallel and diagonal neighbors in 2D or 8 if you have axis-parallel only in 3D).
Quite likely, it is not terribly efficient but neither terribly inefficient to do the slicing manually ( (y[:, :, 1:, 1:] - y[:, :, :-1, :-1]).abs().sum() and friends).
Best regards
Thomas