Following up on my previous email, I have found a workaround to those issues, but now I have a problem with the solvers: Traceback (most recent call last): File "/Users/matak/Downloads/Inequality-constraint-Rockafellar/solvers.py", line 32, in solver_F F_problem = NonlinearVariationalProblem(F_solve, w)#, aP=Ap) File "/Users/matak/firedrake/src/firedrake/firedrake/variational_solver.py", line 47, in __init__ raise ValueError("Provided residual is not a linear form") ValueError: Provided residual is not a linear form There two issues here: 1. I cannot give a preconditioning matrix aP -- is this not possible for nonlinear problems? 2. I suspect the error above is because I have functions that depend on the solution, and I don’t know if that is acceptable in firedrake. I have something like the following: u = Function(V) v = TestFunction(V) maximum = interpolate(0.5*(u + abs(u)), V) F = v*(u + maximum)*dx Best wishes, Anna.