I would recommend trying a change of variables to turn the integration region into a rectangle. For example, let y=u/x, so that the integral becomes
In general, a 2d cubature code can be more efficient than nested 1d quadratures, assuming you are using adaptive quadrature/cubature. This is because the inner 1d integral, if it is performed as an independent 1d adaptive quadrature, can waste a lot of integrand evaluations trying to refine the 1d integral to high relative accuracy even if its overall contribution to the integral is small.
(There are also specialized cubature algorithms for triangular domains, but I’ve never seen an implementation of a nested/adaptive one.)