Over on Slack @avehtari and I have a thread where we recently discussed a path to speed-ups for cases like yours where there is high redundancy in the covariance matrix (in your case thanks to the grid layout of your isotropic space). Code implementing my unique-distances trick is here (though n.b. the proportion-of-variation approach would need tweaking for more appropriate real-world usage), and here’s a plot of the speedups over gp_exp_quad_cov() I get:
Note that the above plot only shows the timings for just the covariance computation, and doesn’t include the more expensive cholesky decomposition time; with that added to both, the max speedup is only about 1.2x. But @avehtari suggests in the Slack thread that the circulant-matrix approach will make the cholesky decomposition (or whatever is used in it’s place; I haven’t clicked through to the examples he provided) should be faster too.
Note that the code at the repo linked above handles the 1D case, but should be easily generalized to your 2D isotropic case.