> On Jun 10, 2021, at 3:29 PM, Kaushik Vijaykumar <[email protected]> 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
