Complexity of a weirdo two-dimensional sorting problem
To elaborate on the comments of Will Sawin and fedja: The question isn't a sorting problem, but it is a matching problem. If $S$ is your arbitrary set and $G = [n]^2$ is your grid, then you are marrying elements of $S$ to elements in $G$, where the happiness of each marriage is your dot product $p \cdot f(p)$. Any happiness function on $S \times G$, not necessarily one that is bilinear in the plane, can be maximized in polynomial time. That's because the convex hull of the set of permutation matrices has few facets: It's the Birkhoff polytope of doubly stochastic matrices. You can apply the general theorem that linear programming with polynomially many facets can be done in polynomial time. For the specific case of the Birkhoff polytope, there is an optimized linear programming algorithm, the Hungarian algorithm, that was discovered and known to be fast before the result that LP is polynomial time in general.
The Hungarian algorithm is already much faster than generic optimization strategies (which generally aren't polynomial time at all). The remaining question is whether you can devise a sorting algorithm that's even faster, maybe even quasilinear time like linear sorting. It's difficult to rule that out, but the setup is sufficiently complicated that I'd be surprised/impressed.
I guess one thing is worth trying: Whether you can prove that you climb to the maximum by separately sorting each row and column of the grid until the permutation is both row-stable and column-stable. That looks a bit like linear programming, and it also looks like a shortcut. Even if this idea doesn't work, it might ably accelerate those linear programming algorithms that travel along the 1-skeleton of the polytope, such as the simplex algorithm. You could also throw in sorting along other straight lines. The Hungarian algorithm does not travel along the 1-skeleton; instead it's a dual/primal algorithm. Actually, traveling along the 1-skeleton of the Birkhoff polytope is a tricky matter because the vertices are highly degenerate. But there are tricks for that such as desingularizing the polytope by making a generic perturbation of the facets.
Here are some more cross-connections: The problem can be massaged into a least-squares matching problem $$\sum_{p\in S} \|(-f(p))-p\|^2 \to \min$$ between $S$ and the grid rotated by $180^\circ$ (which in this case happens to be the same as the original grid, after translation), and hence it becomes a problem of optimal mass transport. The classical reference is:
Franz Aurenhammer, Friedrich Hoffmann, and Boris Aronov, Minkowski-type theorems and least-squares clustering. Algorithmica 20: 61–76 (1998). DOI: 10.1007/PL00009187
There is an implementation in R: https://rdrr.io/cran/transport/man/aha.html
See also my little note Two applications of point matching. I would also be interested to find a provably faster algorithm than the general matching procedures, which would make use of the geometric situation. (Fast practical algorithms and fast approximation algorithms are mentioned in the references of these papers and the software.)