I do not understand how to get asymptotics for the double sum $$\sum_{ab^2 < x} ab^2$$
If I sum over $a$ first, I get $$\sum_{b^2<x} b^2 \sum_{a < x/b^2} a = \frac12 \sum_{b^2<x} b^2 \frac{x^2}{b^4} \sim \frac{\zeta(2)}{2} x^2.$$
If I sum over $b$ first, I get $$\sum_{a<x} a \sum_{b < \sqrt{x/a}} b^2 = \frac13 \sum_{a<x} a \left( \frac{x}{a} \right)^{3/2} \sim \frac{2}{3} x^2. $$
What is happening. Am I doing something forbidden in switching summations this way? I suppose it is a possible divergence problem, so I do I get precise asymptotics for this sum?
TLDR: Your second method uses two asymptotics. The asymptotic approximation for the inner summation throws out "lower order terms" which end up mattering for outer summation.
For integers $N \ge 0$, we have $$\displaystyle\sum_{k = 1}^{N}k^2 = \dfrac{N(N+1)(2N+1)}{6}.$$ Combining this with the inequality $y-1 \le \lceil y \rceil-1 \le y$, we get $$\dfrac{(y-1)y(2y-1)}{6} \le \displaystyle\sum_{k = 1}^{\lceil y \rceil-1}k^2 \le \dfrac{y(y+1)(2y+1)}{6},$$ i.e. $$\dfrac{1}{3}y^3-\dfrac{1}{2}y^2+\dfrac{1}{6}y \le \displaystyle\sum_{k < y}k^2 \le \dfrac{1}{3}y^3+\dfrac{1}{2}y^2+\dfrac{1}{6}y$$ for any real number $y \ge 1$. Applying this yields $$\dfrac{x^{3/2}}{3a^{3/2}}-\dfrac{x}{2a}+\dfrac{x^{1/2}}{6a^{1/2}} \le \sum_{b < \sqrt{x/a}}b^2 \le \dfrac{x^{3/2}}{3a^{3/2}}+\dfrac{x}{2a}+\dfrac{x^{1/2}}{6a^{1/2}},$$ and thus, $$\sum_{a < x}\left[\dfrac{x^{3/2}}{3a^{1/2}}-\dfrac{x}{2}+\dfrac{x^{1/2}a^{1/2}}{6}\right] \le \sum_{a < x}a\sum_{b < \sqrt{x/a}}b^2 \le \sum_{a < x}\left[\dfrac{x^{3/2}}{3a^{1/2}}+\dfrac{x}{2}+\dfrac{x^{1/2}a^{1/2}}{6}\right].$$
However, note that $\displaystyle\sum_{a < x}\dfrac{x^{3/2}}{3a^{1/2}} \sim \dfrac{2}{3}x^2$, $\displaystyle\sum_{a < x}\dfrac{x}{2} \sim \dfrac{1}{2}x^2$, and $\displaystyle\sum_{a < x}\dfrac{x^{1/2}a^{1/2}}{6} \sim \dfrac{1}{9}x^2$. So the lower order terms that the approximation $\displaystyle\sum_{k < y}k^2 \sim \dfrac{1}{3}y^3$ drops end up contributing a significant amount to the double sum. Hence, all you know from the second method is that the double sum is asymptotically bounded between $\dfrac{5}{18}x^2$ and $\dfrac{23}{18}x^2$.
We can analyze the first method in the same manner. Using $$\dfrac{1}{2}y^2-\dfrac{1}{2}y \le \sum_{k < y}k \le \dfrac{1}{2}y^2+\dfrac{1}{2}y,$$ we get $$\dfrac{x^2}{2b^4}-\dfrac{x}{2b^2} \le \sum_{a < x/b^2}a \le \dfrac{x^2}{2b^4}+\dfrac{x}{2b^2},$$ and thus, $$\sum_{b < \sqrt{x}}\left[\dfrac{x^2}{2b^2}-\dfrac{x}{2}\right] \le \sum_{b < \sqrt{x}}b^2\sum_{a < x/b^2}a \le \sum_{b < \sqrt{x}}\left[\dfrac{x^2}{2b^2}+\dfrac{x}{2}\right].$$ As you already noted, $\displaystyle\sum_{b < \sqrt{x}}\dfrac{x^2}{2b^2} \sim \dfrac{\zeta(2)}{2}x^2$. The other term scales like $\displaystyle\sum_{b < \sqrt{x}}\dfrac{x}{2} \sim \dfrac{1}{2}x^{3/2} = o(x^2)$. Hence, the double sum does indeed scale like $\dfrac{\zeta(2)}{2}x^2$.