I would suggest creating your own solver. PETScSNESSolver allows you to provide the lower and upper bounds for your function.
So the way forward could be to create a subclass of NonlinearProblem, calculate the jacobian as you have above, and then use PETScSNESSolver to solve the problem like here, along with providing the bound vectors. Basically the array of dofs for your lower and upper functions (namely lower.vector() and upper.vector()). See also here
Edit: It does seem that you can do the same with NonlinearVariationalProblem, I can take a look at it sometime later today.