Non linear minimisation on constrained boundary

Hi all,

I have a question about bounded non linear PDEs with a mixed element state, in this case of a functional including a \phi(\dot{p}) plastic dissipation term, shown below. Instead of using a ufl.conditional with \phi(\dot{p}<0) = +\infty to enforce the bounds, I’ve used SetVariableBounds with a non linear SNES solver.
image

The problem begins with \dot{p} = 0 , i.e. on the edge of the boundary, hence its derivative is undefined so I believe that purely gradient based solvers may not converge at the boundary, as is the case with my current setup. Is this the case and are there any solvers that handle problems at boundaries well?

Kind regards,

Magnus

These are the solver conditions, for reference.

snes = PETSc.SNES().create()
opts = PETSc.Options()
opts['snes_converged_reason'] = None
#opts['snes_linesearch_type'] = 'basic'
opts['snes_monitor'] = None
#opts['snes_linesearch_monitor'] = None
snes.setFromOptions()
snes.setType('vinewtonrsls')
snes.setFunction(problem.F, b)
snes.setJacobian(problem.J, J)
snes.setTolerances(atol=1e-10,rtol=1.0e-7, max_it=20)
snes.getKSP().setType("preonly")
snes.getKSP().setTolerances(rtol=1.0e-7)
snes.getKSP().getPC().setType("lu")
snes.getKSP().getPC().setFactorSolverType("mumps")