This was a pain to debug. The fix is the following:
w = df.Function(W)
pv, sv = df.TestFunctions(W)
phi, psi = df.split(w)
In debugging I noticed that the matrix you’re assembling is empty. This implies that the derivative isn’t correctly formed. And indeed this lead to finding that w was defined as a Function and not actually used in your nonlinear formulation. Therefore the Gateaux derviative is zero, and no matrix is assembled.
The fix will allow you to assemble the problem; however, with regards to getting the solver to work I refer you to some tips and tricks. The nonlinear problem you’ve written is not trivial given the fractional powers.