[petsc-users] Fwd: Implementing a homotopy solver
zakaryah
zakaryah at gmail.com
Sun Jul 22 20:09:57 CDT 2018
Thanks Matt and Barry,
Matt - if I do the calculation in FormJacobian(), which makes by far the
most sense and is as per your suggestion, do I need to set the operators of
the SNES's KSP back to whatever they were before I set them? The linear
system I want to solve within FormJacobian() involves the Jacobian matrix
itself, and I want to remove the "nullspace" from that same matrix within
FormFunction().
Barry - I'm trying to implement a homotopy solver. In short, I have a
system of n nonlinear equations in n variables, F(x), which is hard to
solve because the Jacobian tends to become singular using standard
methods. I want to add an auxiliary variable, lambda, to create a
homotopy: H(lambda,x) = lambda*F(x) + (1-lambda)G(x), where G is "easy to
solve", and the idea is that the Jacobian of the n+1 variable system will
not become singular along the curve H(lambda,x) = 0.
The method involves adding an equation to H so that the Jacobian H' is
square. The "submatrix" refers to the n x (n+1) matrix which represents
the Jacobian without the added equation, whereas my FormJacobian() forms
the entire (n+1) x (n+1) matrix H'. I only refer to the submatrix because
it has a nullspace, and I want to find it by solving a linear system
designed for this purpose, H' u = b, where b is not the zero vector. H'
has no nullspace, but I want to remove the projection of u from my SNES
solution vector, as u IS in the nullspace of the submatrix.
The RHS vector b can be calculated outside the SNES solve. I guess my
FormJacobian() should look like this:
FormJacobian(SNES snes, Vec x, Mat Amat, Mat Pmat, void *ctx) {
KSP ksp;
MatNullSpace unull;
user_struct *user = (user_struct*)ctx;
calculate Amat
assemble Amat
if (Pmat != Amat) {
assemble Amat
}
SNESGetKSP(snes,&ksp);
KSPSetOperators(ksp,Amat,Pmat);
KSPSolve(ksp,user->b,user->u);
MatNullSpaceCreate(PetscObjectComm((PetscObject)Amat), PETSC_FALSE, 1,
&(user->u),&unull);
MatSetNullSpace(Amat,unull);
MatNullSpaceDestroy(&unull);
}
Does this look right?
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180722/3ddfd2fa/attachment.html>
More information about the petsc-users
mailing list