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
_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake