Counting Total Number of Non-Equivalent Configurations in a 2-D Grid
We will use Burnside's lemma, plus some additional optimizations.

Note that row swaps and column swaps commute, so the group of allowable transformations to the grid is exactly $G = S_w \times S_h$, where $S_n$ is the symmetric group on $n$ symbols.  Note that $G$ embeds in $S_{wh}$ via a canonical embedding $\varphi : G \to S_{wh}$, given by $\varphi(\pi,\sigma) = \tau$ where $\tau$ is defined by $\tau(i,j) = (\pi(i), \sigma(j))$.  Therefore, we will think of $G$ simultaneously as the direct product $S_w \times S_h$ and also as a subgroup of $S_{wh}$.

Now let's count the number of grids that are fixed by a group element $g=(\pi,\sigma) \in G$.  Define the cycle type of a permutation $\pi$ to be the multiset of lengths of the cycles in the cycle decomposition of $\pi$.  If $\pi \in S_n$, the cycle type of $\pi$ is a partition of $n$.  With this in mind, we will let $\ell(\pi)$ denote the cycle type of $\pi \in S_w$ and $\ell(\sigma)$ the cycle type of $\sigma \in S_h$.  Finally, let $\ell((\pi,\sigma))$ denote the cycle type of $(\pi,\sigma)$, thinking of $(\pi,\sigma)$ as an element of $S_{wh}$.  We can express $\ell((\pi,\sigma))$ in terms of $\ell(\pi)$ and $\ell(\sigma)$: let $x,y$ range over $x \in \ell(\pi), y \in \ell(\sigma)$, and for each $x,y$, $\ell((\pi,\sigma))$ has $\gcd(x,y)$ copies of $\operatorname{lcm}(x,y)$.

Now the number of grids that are fixed by $(\pi,\sigma)$ is

$$f(\pi,\sigma) = k^{|\ell(\pi,\sigma)|},$$

where $|\ell(\pi,\sigma)|$ counts the number of cycles in $(\pi,\sigma) \in G$ (multiplicities are counted).  By the above characterization of $\ell(\pi,\sigma)$, we see that

$$f(\pi,\sigma) = k^{\sum_{x \in \ell(\pi)} \sum_{y \in \ell(\sigma)} \gcd(x,y)}.$$

To apply Burnside's lemma, it remains to sum $f(g)$ over all elements $g \in G$, as Burnside's lemma tells us that the total number of non-equivalent grids is

$$N = {1 \over |G|} \sum_{g \in G} f(g) = {1 \over w! h!} \sum_{g \in G} f(g).$$

In principle, we could compute this sum by enumerating all $w! h!$ elements in $G$, then evaluating $f$ on each, and summing the result.

It turns out this can be optimized further, as $f((\pi,\sigma))$ depends only on $\ell(\pi)$ and $\ell(\sigma)$, but not on any other aspect of $\pi$ or $\sigma$.  Thus, we can more efficiently evaluate the sum using the formula

$$\sum_{g \in G} f(g) = \sum_{P,Q} n(P,Q) \times k^{\sum_{x \in P} \sum_{y \in Q} \gcd(x,y)},$$

where the inner sum ranges over all partitions $P$ of $w$ and all partitions $Q$ of $h$, and where we define

$$n(P,Q) = |\{(\pi,\sigma) \in G : \ell(\pi)=P, \ell(\sigma)=Q\}|$$

to count the number of group elements of $G$ whose cycle types are $P,Q$.  Now we have

$$n(P,Q) = n_w(P) \times n_h(Q),$$

where $n_w(P) = |\{\pi \in S_w : \ell(\pi)=P\}|$ is the number of permutations in $S_w$ with cycle type $P$, and similarly $n_h(Q)$ is the number of permutations in $S_h$ with cycle type $Q$.  Finally, it is easy to compute $n_w(P)$ for any desired partition $P$; it is

$$n_w(P) = {w! \over \prod_{x \in P} x}.$$

Plugging in, we find that the total number of non-equivalent grids is given by

$$N = \sum_{P,Q} {k^{\sum_{x \in P} \sum_{y \in Q} \gcd(x,y)} \over \prod_{x \in P} x \prod_{y \in Q} y}.$$

This can be computed by two nested for-loops that iterate over all partitions $P$ of $w$ and all partitions $Q$ of $h$.

What is the running time of this algorithm?  The  asymptotic running time is approximately $O(p(w) p(h))$, where $p(n)$ is the number of partitions of the integer $n$.  It is known that asymptotically $p(n) = O(\exp\{\pi \sqrt{2n/3} \})$, so it follows that the asymptotic running time of this algorithm is approximately $O(\exp\{\pi \sqrt{2/3} (\sqrt{w}+\sqrt{h})\})$, or equivalently, approximately $O(2^{3.7(\sqrt{w}+\sqrt{h})})$.  This should allow you to handle the full range of values $1 \le w,h \le 12$, as for that range we have $\sqrt{w}+\sqrt{h} &lt; 7$ and $2^{3.7 \times 7} \approx 2^{26}$ is not too large.