Hi Justin, On 27/06/16 23:35, Justin Chang wrote:
Hi all,
So I am attempting to use TAO's mixed complementarity solvers and it basically requires a formGradient/Constraints routine and a formJacobian routine.
I am solving a variational inequality via TAO and want the linear operator A to be my jacobian. However, I don't think it was correctly done because my semi smooth linear solver does not converge (but it should!).
Attached is the code. Line 108 is what TaoSetJacobianRoutine() calls. Basically I want to assemble the linear operator into the matrix petsc_J. Is what I have the correct implementation?
...
# Complementarity solver def formJac(tao, petsc_x, petsc_J, petsc_JP, J=None, b=None, a=None, bcs=None):
Here is your problem. You need to be filling the entries of petsc_J and maybe petsc_JP. Instead, you are reassigning to the variable name, but not modifying the entries. If this is like SNES, where PETSc just expects you to assemble into the Jacobian you passed in using setJacobian, then I think you want to do: def formJac(tao, petsc_x, petsc_J, petsc_JP, J=None, b=None, a=None, bcs=None): # Reassemble into the tensor J. J = assemble(a, bcs=bcs, tensor=J) J.M._force_evaluation() # Check that we're getting it right. assert petsc_J == J.M.handle ...
taoSolver.setJacobian(formJac,J.M.handle, kargs={'J': J, 'b': b, 'a': a, 'bcs': bcs})
Use J._M.handle here (since you don't need to evaluate it now). Then I think this should now work. Cheers, Lawrence