[petsc-users] Copy stiffness matrix to jacobian matrix in form jacobian

[Kaushik Vijaykumar]

Thanks Barry for your comments. Based on your suggestion I tried the
following

Main()
{
ierr = SNESSetFunction(snes,r1,FormFunction,&this_ctx);CHKERRQ(ierr);
ierr = SNESSetJacobian(snes,jac,NULL,FormJacobian,&this_ctx);CHKERRQ(ierr);
ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
ierr = SNESSolve(snes, NULL, u);CHKERRQ(ierr);
}

In formfunction I populate the jacobian "*jac*"

Formfunction(snes,Vec x,Vec F,void* ctx)
{
ierr = SNESGetJacobian(snes,&jac,Null,NULL,NULL);CHKERRQ(ierr);

// "*Populate jac*"
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

}
FormJacobian(SNES snes,Vec x,Mat jac,Mat B,void *ctx)
{
ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);

ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}

When I execute this I get a Petsc error, object in wrong state missing
diagonal entry.

I am not sure what I am missing here?

Thanks
Kaushik


On Thu, Jun 10, 2021 at 6:37 PM Barry Smith <bsmith at petsc.dev> wrote:

>
>
> > On Jun 10, 2021, at 3:29 PM, Kaushik Vijaykumar <kaushikv318 at gmail.com>
> wrote:
> >
> > Hi everyone,
> >
> > I am trying to copy the stiffness matrix that I generated in form
> function to the jacobian matrix jac in form jacobian. The piece of code for
> that:
> >
> >
>    If you have created jac and B in your main program and passed them with
> SNESSetJacobian() then you should be able to just use
>
> > ierr = MatCopy(ptr->K,jac,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
>
> Generally jac and B are the same matrix so you don't need the second copy.
> > ierr = MatCopy(ptr->K,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
>
> If they are sometimes the same and sometime not you can do
>
> if (jac != B) {ierr =
> MatCopy(ptr->K,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);  }
>
>
>    The following won't work because you are creating new local matrices in
> your form jacobian that don't exist outside of that routine.
>
> >
> > ierr = MatDuplicate(ptr->K,MAT_COPY_VALUES,&jac);CHKERRQ(ierr);
> > ierr = MatDuplicate(ptr->K,MAT_COPY_VALUES,&B);CHKERRQ(ierr);
> >
> > ierr = MatCopy(ptr->K,jac,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
> > ierr = MatCopy(ptr->K,B,SAME_NONZERO_PATTERN);CHKERRQ(ierr);
> >
> >
>    If the matrix you are creating in form function is the Jacobian then
> you don't need to use the memory and time to make copies. You can just form
> the matrix directly into the Jacobian you will use for SNES. You can obtain
> this with a call
>
>    SNESGetJacobian(snes,&jac,NULL,NULL,NULL);
>
> in your form function and just put the matrix values directly into jac.
> Then your form Jacobian does not need to do anything but return since the
> jac already contains the Jacobian.
>
> On the off-change you are using the iteration A(x^n) x^{n+1} = b(x^n} and
> not using Newtons' method (which is why you compute A in form function,
> then you might be better off using SNESSetPicard() instead of
> SNESSetFunction(), SNESSetJacobian().
>
> >
> > Thanks
> > Kaushik
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210611/5f8cf792/attachment-0001.html>