[petsc-users] Copy stiffness matrix to jacobian matrix in form jacobian
Barry Smith
bsmith at petsc.dev
Mon Jun 14 14:45:43 CDT 2021
It is best to cut and paste the entire error message, not just two lines. We are missing all context to see where the error is occurring and possible causes.
> On Jun 14, 2021, at 6:20 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"
>
>
> Thanks
> Kaushik
>
> On Fri, Jun 11, 2021 at 11:33 AM Barry Smith <bsmith at petsc.dev <mailto: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 <mailto: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 <mailto:bsmith at petsc.dev>> wrote:
>>
>>
>> > On Jun 10, 2021, at 3:29 PM, Kaushik Vijaykumar <kaushikv318 at gmail.com <mailto: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/7ad41a2e/attachment-0001.html>
More information about the petsc-users
mailing list