NewtonSolver update solution with inequality constraint

Question

Does anyone have an idea to update the solution in Newton iteration scheme to constrain the solution to be non-negative. I know you can do something similaer with the PETScSNESSolver but I'd like to implement my own constraint algorithm. I have tried something like

class MyNonNegativeSolutionNewtonSolver(NewtonSolver):

def __init__(self, u):
       NewtonSolver.__init__(self)

def update_solution(self, x, dx):
    i = (x-dx)<0.0
    dx[i] = 0.0
    NewtonSolver.update_solution(self,x,dx)

But I am concern that this might lead to some oscillations as the solution approaches zero probably because my Jacobian is not constrained accordingly. Any suggestions on how to implement solutions constrained with inequalities are welcome.

Answer 1

You cannot do a simple projection to impose the constraint. With the constraint the problem to solve changes drastically. I strongly suggest you to use the specifically designed algorithm available in PETScSNES and PETScTAO. You can use them in python either through the dolfin interface or via petsc4py