After several workarounds, I finally found a solution that works nicely!
The general issue was that, in my variational form, there are terms of the form
\ln\left(\frac{u}{1-u}\right).
The Newton solver sometimes produced negative values of u as a numerical artifact, which led to the failure of the solver.
I got rid of diverging solutions by forcing the argument of the logarithm to remain positive (DOLFIN_EPS) with a conditional statement for u>0:
ln ( conditional( gt(u, 0), u/(1-u), DOLFIN_EPS ) )
I think that this is rather straightforward, but I didn’t think of it. I had the idea when I saw this post:
So thank you, @kamensky!