I think this can be answered by drawing from sequential Monte-Carlo, and quasi Monte-Carlo methods.
The key concept that underlies your problem is the exploration-exploitation dilemma:
- Exploration: you want to cover as much as possible of the space where $f$ is defined;
- Exploitation: but ideally you would rather spend time sampling places where the value of $f$ varies a lot, rather than where it is flat.
For the exploration part, quasi Monte-Carlo methods tell us to choose low-discrepancy sequences in order to cover the space as efficiently as possible. In 2D the method you choose does not matter too much; you basically need a grid-like structure, see the Sobol sequence for example. If you just want to explore as much as possible, this answers your question (see the Koksma-Hlawka inequality).
For the exploitation part, we need to adopt a more dynamic perspective; somehow, we need to learn from the samples that already exist, in order to focus on places that matter to us, as we sample more and more. Here, importance sampling gives us a way of focusing on places of interest. For example if we were interested in finding the peaks of $f$, given a set of $N$ existing samples $(x_i, f(x_i))$, we would resample in the neighbourhood of $x_k$ with a probability proportional to $f(x_k)$. In your case, you may want to resample where the gradient of $f$ is large for example.
How to balance exploration and exploitation is a question without an answer; it really depends on your problem and on the resources available.
To give a concrete answer to your problem, you could consider a 3x3 grid covering your 2D space, with a sample at the centre of each cell. From that, you could estimate the "gradient" in each cell, by comparing the value obtained with the neighbours. Iteratively then:
- select a subset of cells with large gradients (and possibly other cells selected at random),
- divide them further with a nested 3x3 grid;
- evaluate $f$ at the centre of each smaller cell obtained;
- re-evaluate the "gradient" in each cell by comparing the value obtained with the neighbours.
This would be an adaptive way of sampling $f$. Obviously this is just a sketch of a solution, and many details need to be addressed (mainly how to compare values obtained between boxes of different sizes, and what is the variance to be expected on the estimate of the gradient).