It's fun to be square

With some easy integration and some fun with the polygamma function: $\begin{array}{rcl} S \; = \; \displaystyle \int_0^1 \int_0^1 \frac{y}{1-x^3y^3}\,dy\,dx & = & \displaystyle \int_0^1 \left(\int_0^1 \frac{y}{1 - x^3y^3}\,dx\right) \,dy \\ & = & \displaystyle \int_0^1 \left(\int_0^y \frac{1}{1-u^3}\,du\right)\,dy \; = \; \int_0^1 \left(\int_0^y \sum_{n=0}^\infty u^{3n}\,du \right)\,dy \\ & = & \displaystyle \int_0^1 \sum_{n=0}^\infty \frac{y^{3n+1}}{3n+1}\,dy \; = \; \sum_{n=0}^\infty \frac{1}{(3n+1)(3n+2)} \\ & = & \displaystyle \sum_{n=0}^\infty \left(\frac{1}{3n+1} - \frac{1}{3n+2}\right) \; = \; \tfrac13\sum_{n=0}^\infty \left(\frac{1}{n + \frac13} - \frac{1}{n+\frac23}\right) \\ & = & \displaystyle \lim_{n \to \infty} \frac13\Big[\psi\big(n +\tfrac13\big) - \psi\big(\tfrac13\big) - \psi\big(n + \tfrac23\big) + \psi\big(\tfrac23\big)\Big] \\ & =& \displaystyle \tfrac13\big[\psi\big(\tfrac23\big) - \psi\big(\tfrac13\big)\big] \; = \; \frac{\pi}{3\sqrt{3}} \end{array}$ making the answer $(3\sqrt{3})^2 = \boxed{27}$ .

As the limits of integration are constant we can exchange $dx$ and $dy$ .

So the integral becomes:-

$\Large \int_{0}^{1}\int_{0}^{1} \frac{y}{1-(xy)^{3}} dx dy$

Substituting $xy=t$

we have:-

$\Large \int_{0}^{1}\int_{0}^{y} \frac{1}{1-t^{3}} dt dy$

Now we want to change the order of $dt$ and $dy$

So we have:-

$\Large \int_{0}^{1}\int_{t}^{1} \frac{1}{1-t^{3}} dy dt$

Evaluating we have:-

$\Large \int_{0}^{1}\ \frac{1-t}{1-t^{3}}dt$

$\Large = \int_{0}^{1}\frac{1}{1+t+t^{2}}dt$

$\Large = \int_{0}^{1}\frac{1}{(t+\frac{1}{2})^{2} + \frac{3}{4}}dt$

which is

$\frac{2}{\sqrt{3}}(\tan^{-1}(\sqrt{3}) - \tan^{-1}(\frac{1}{\sqrt{3}}))$

$=\frac{\pi}{3\sqrt{3}}$

Also we could do it using series:-

The series that @Mark Hennings arrived at can be solved using complex numbers too and using taylor series of $ln(1-x)$

Using that it can be written as

$\Re(\frac{1}{3}(\omega^{2}-\omega)(ln(1-\omega)-ln(1-\omega^{2})))$

where omega is a complex cube root of unity and $\Re(.)$ denotes the real part.

I'll leave it upto the really interested to find out how the expression in complex numbers written above evaluates to the series

$\sum_{r=0}^{\infty}(\frac{1}{3r+1} - \frac{1}{3r+2})$