This depends on exactly what you mean by non-uniform, or (equivalently) on how big the loop is. In particular, the important criterion is whether the field changes appreciably over distances that are about the same size as the loop.
If the field changes throughout space, but the loop is small enough that the field doesn't change much from point to point on the loop, then the uniform-field formula $\boldsymbol\tau=\boldsymbol\mu\times\mathbf B(\mathbf r)$ still applies. In essence, the field is locally uniform, though the direction and magnitude it's uniform on can change from place to place.
If the loop is big enough that the field changes appreciably over its span then there's nothing for it but to integrate the local torque on each bit of circuit and add them up, which gives you
$$\boldsymbol\tau
=\oint_\mathcal{C}\mathbf r\times\mathbf F(\mathbf r)\:\mathrm d l
=\oint_\mathcal{C}\mathbf r\times(\hat{\mathbf t}I\times\mathbf B(\mathbf r))\:\mathrm d l
=\oint_\mathcal{C}\mathbf r\times(I\mathrm d \mathbf l\times\mathbf B(\mathbf r)).
$$
There really isn't much you can do to simplify it beyond that without special assumptions. The integral is a line integral, of exactly the same sort you use to calculate the magnetic dipole moment $\boldsymbol \mu$ itself.
The following is a bit more technical and uses a standard amount of vector calculus as used in undergraduate electromagnetism courses.
If the field were constant, then you can manipulate it more freely to get
\begin{align}
\boldsymbol\tau
=\oint_\mathcal{C}\mathbf r\times(I\mathrm d \mathbf l\times\mathbf B)
=I\oint_\mathcal{C}\left[(\mathbf B·\mathbf r)\mathrm d\mathbf l-(\mathbf r·\mathrm d\mathbf l)\mathbf B\right],
\end{align}
using a vector triple product. Here the second integral vanishes, because Stokes' theorem implies that
$$
\oint_\mathcal{C}\mathbf r·\mathrm d\mathbf l=\int_\mathcal{S}(\nabla\times\mathbf r)·\mathrm d\mathbf S=0.
$$
You're left with the first integral, $\boldsymbol\tau=I\oint_\mathcal{C}(\mathbf B·\mathbf r)\mathrm d\mathbf l$, which in component notation reads $\tau_i=IB_j\oint_\mathcal{C} x_j \mathrm dx_i$ (I use Einstein summations throughout). As it turns out, the indices on that last integral are antisymmetric, i.e.
$$\oint_\mathcal{C} x_j \mathrm dx_i=-\oint x_i \mathrm dx_j,$$
which you can again prove via Stokes' theorem since
\begin{align}
\oint_\mathcal{C}\left[ x_j \mathrm dx_i+x_i\mathrm dx_j\right]
=\oint_\mathcal{C}\left[ x_j (\nabla x_i)·\mathrm d\mathbf l+x_i(\nabla x_j)·\mathrm d\mathbf l\right]
=\oint_\mathcal{C}\nabla(x_ix_j)·\mathrm d\mathbf l
=\int_\mathcal{S}\nabla\times\nabla(x_ix_j)·\mathrm d\mathbf S
=0.
\end{align}
What this means is that you can "fold" the old integral into two antisymmetric parts, as
$$\tau_i=IB_j\oint_\mathcal{C} \frac{x_j \mathrm dx_i-x_i \mathrm dx_j}{2},$$
and the antisymmetrized integral is now exactly the magnetic dipole moment
$$
\boldsymbol\mu=\frac I2\oint_\mathcal C\mathbf r\times\mathrm d\mathbf l
$$
which in component notation reads
$$
\mu_k=I\oint_\mathcal C\varepsilon_{kij}x_i\mathrm dx_j
,\quad\text{or in other words}\quad
\mu_k\varepsilon_{kij}=I\oint_\mathcal C \left[x_i\mathrm dx_j-x_j\mathrm dx_i\right].
$$
Back to the torque this means
$$\tau_i=\frac12 B_j\varepsilon_{kji}\mu_k,$$
or in vector notation
$$
\boldsymbol\tau=\boldsymbol\mu\times\mathbf B.
$$
So what is it I've done? All of this work has been to take an integral that had the magnetic field inside it, and factorize it into a system-dependent part and a field dependent part:
$$
\boldsymbol\tau
=\oint_\mathcal{C}\mathbf r\times(I\mathrm d \mathbf l\times\mathbf B)
=\left(\frac I2\oint_\mathcal C\mathbf r\times\mathrm d\mathbf l\right)\times\mathbf B.
$$
That's really what the magnetic dipole moment $\boldsymbol\mu$ really is.
OK, sorry, that got out of hand, but I'll get back on topic now. What does this have to do with nonuniform magnetic fields? Well, if your field varies very slowly with respect to the dimensions of the loop, you can suppose that $\mathbf B$ is constant in your integral, and you get the calculation above. The next thing you might try is suppose that $\mathbf B$ is almost constant throughout the extent of the loop, but that you do need to consider the first-order term of its Taylor series. Thus, you could suppose that
$$
\mathbf B(\mathbf r)=\mathbf B(\mathbf r_0)+(\mathbf r·\nabla)\mathbf B(\mathbf r_0),
$$
or more clearly in component notation
$$
B_i(\mathbf r)=B_i(\mathbf r_0)+x_k\frac{\partial B_i}{\partial x_k}(\mathbf r_0).
$$
You can then do the same game I've done above with the second term, with the added complication that you have an extra factor of $x_k$ in your integral. The result will be something which depends on the first-order derivative of the field, $\frac{\partial B_i}{\partial x_k}(\mathbf r_0)$, and if you're clever you can factorize out the dependence on the loop into a single factor. This factor will in general be (i) a matrix, or in fancy-speak a tensor, (ii) the integral of a homogeneous second-degree polynomial over the loop, and (iii) it is called the quadrupole moment of the loop.
In general, such calculations are long and messy, but if you're feeling brave I encourage you to have a go at deriving an expression for the quadrupole moment and seeing how simple of an expression you can get for the moment itself and its interaction with the magnetic field's gradient to get the torque. As a start, though, here's one expression for the quadrupolar part of the torque, which I encourage you to derive
$$
\tau_i=\frac{\partial B_m}{\partial x_k}·I\oint_\mathcal{C}\left[
x_kx_m\mathrm d x_i-\delta_{im}x_k x_j\mathrm dx_j
\right].
$$
Of course, there's only so much you can do with only first-order derivatives. If your field varies slightly faster than that - or if your circuit is somewhat bigger - then the next thing you can try is a second-order Taylor expansion of the field. This gives you a third term which depends on the second derivatives of the magnetic field and on a system-dependent term which is (i) an unwieldy object with three different indices, called a rank-3 tensor, (ii) the integral of a homogeneous third-degree polynomial over the loop, and (iii) is called the octupole moment of the loop. And after that, you can go to even higher orders, and on and on it goes until you feel your calculation is accurate enough, or you give up through exhaustion.
I can tell you, though - none of the formulas there are pretty. That's why they're hard to find online.
Parallel currents attract. The solenoid, with current in the direction shown, has current parallel to that in the loop, and the loop will be attracted to the solenoid. A uniform B field, however, would put no net force on such a loop; any force on one section of the loop would be cancelled by equal and opposite force on
another section of the loop (the opposite side of the loop carrying the same
current in the opposite direction).
The radial (radiating as if it were coming from a point in the center) B field
component is not uniform (it's strong in the center, and points in various
directions). That part of the B field puts a net force on the current-carrying ring.
Best Answer
How to do it the hard way
The Lorentz-force on a curve $C$ carrying a uniform current I is given by the line integral: $$\vec F = I \int_C d\vec \ell \times \vec B$$ For a constant $\vec B$-field over a closed loop, this reduces to: $$ \vec F = I \left[\oint_C d\vec r \right]\times \vec B = I \vec 0 \times \vec B = \vec 0.$$ The key factor here is that the line integral of $1$ over a closed loop is $\vec 0$ no matter the shape of the loop. It is easiest to see this component-by-component; the $x$-component for example is $\oint_C dx = \Delta x,$ but any net displacement in the positive direction needs to eventually be compensated by net displacement in the negative direction, to bring it back to where it started.
Magnetic force on a current-carrying wire
Why the heck is that first expression true? You just have to be really meticulous about what you mean, basically. A general current loop is a parametric curve $C = \{~\vec r(s) \text{ for all } s \text { in } (0, S)~\}.$ Between two points $s, s+ds$ you'll see a tangent vector $d\vec r = \vec r'(s) ~ds.$ There is probably a linear charge density for the moving charges $\lambda(s)$ over the loop, and they move at a speed $v(s)$ such that $\lambda(s)~v(s) = I(s)$ is the current through the wire: we usually want this to be a constant, so charge isn't "building up" at any point, otherwise that creates $\vec E$ fields which oppose the $\vec B$-field forces involved.
We don't really have a time coordinate as we're only looking at one moment in time, but let me take a typical charge flow over parameter-change $ds$ as happening over a real-time-change $dt$. Then, identifying $v(s)$ as $\left|\frac {d\vec r}{dt}\right| = \left|\vec r'(s)\right| \cdot \left|\frac{ds}{dt}\right|$ we get $\frac{ds}{dt} = I/\big[\lambda(s)~ |\vec r'(s)|\big].$ So that gives an explicit notion for $dt$.
The net Lorentz force due to the magnetic field is of course $$\vec F = \oint_C dq~\vec v\times\vec B, $$and we identify $dq = \lambda(s) ~ |\vec r'(s)|~ ds$ and $\vec v = \frac {d\vec r}{dt} = \vec r'(s) ~\frac{ds}{dt}$ to turn this into: $$\vec F = \int_0^S ds~I(s) ~ \vec r'(s) \times\vec B(s) = I \oint_C d\vec r\times \vec B.$$ So you just get a straightforward line integral.
Another simple way to see this is as the curvy generalization of the well-known result for a length $L$ of wire oriented in direction $\hat n$ with current $I$ going in the $\hat n$ direction: then the force on the wire is $\vec F = L~I~\hat n\times\vec B.$ (Some people write $\vec I = \hat n ~I.$)