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

[Matthew Knepley]

On Mon, Jun 14, 2021 at 7:21 AM Kaushik Vijaykumar <kaushikv318 at gmail.com>
wrote:

> 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"
>

It must be that the matrix you are setting the diagonal for is not the same
matrix as referenced in the error message.

  Thanks,

     Matt


> 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
>>>
>>>
>>

-- 
What most experimenters take for granted before they begin their
experiments is infinitely more interesting than any results to which their
experiments lead.
-- Norbert Wiener

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210614/8c4a2a94/attachment.html>