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

[Kaushik Vijaykumar]

Thanks for the suggestion Barry, I did do that by setting it in multiple
ways first by setting

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

.....


ierr = MatGetOwnershipRange(jac,&istart,&iend);

for(i=istart;i<iend;i++)
{
val = 0.0;
ierr = MatSetValue(jac, i, i, val, ADD_VALUES);CHKERRQ(ierr);
}
ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);
}
I get the same error as before for a single processor run.
"[0]PETSC ERROR: --------------------- Error Message
--------------------------------------------------------------
[0]PETSC ERROR: Object is in wrong state
[0]PETSC ERROR: Matrix is missing diagonal entry 0"


Thanks
Kaushik

On Fri, Jun 11, 2021 at 11:33 AM Barry Smith <bsmith at petsc.dev> wrote:

>
> When I execute this I get a Petsc error, object in wrong state missing
> diagonal entry.
>
>
>    When do you get this error? During a TSSolve? While KSP is building a
> preconditioner.
>
>    Some parts of PETSc require that you put entries (even if they are
> zero) on all diagonal entries in the Jacobian. So likely you need to just
> make sure that your populate jac routine puts/adds 0 to all diagonal
> entries.
>
>   Barry
>
>
>
> On Jun 11, 2021, at 10:16 AM, Kaushik Vijaykumar <kaushikv318 at gmail.com>
> wrote:
>
> 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/20210614/4b7db989/attachment.html>