[petsc-users] Changing the jacobian matrix during SNES iterations

[David Knezevic]

On Wed, Jul 18, 2018 at 3:34 PM, Matthew Knepley <knepley at gmail.com> wrote:

> On Wed, Jul 18, 2018 at 3:25 PM David Knezevic <david.knezevic at akselos.com>
> wrote:
>
>> On Wed, Jul 18, 2018 at 1:59 PM, Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>>> On Wed, Jul 18, 2018 at 1:31 PM David Knezevic <
>>> david.knezevic at akselos.com> wrote:
>>>
>>>> I'm using SNES for a finite element contact solve, in which the
>>>> sparsity pattern of the jacobian can change from one Newton iteration to
>>>> the next (since the nodes on the contact surface move).
>>>>
>>>> In order to handle this I figured the best way would be to destroy the
>>>> jacobian matrix and re-allocate it with a new sparsity pattern inside each
>>>> call to FormJacobian, does that seem like a reasonable approach in this
>>>> context?
>>>>
>>>
>>> Yes.
>>>
>>>
>>>> Also, I recall from an earlier discussion that this matrix
>>>> re-allocation inside FormJacobian is supported by SNES, but I just wanted
>>>> to confirm that?
>>>>
>>>
>>> Yes.
>>>
>>>
>>>> Also, I was wondering if there is any example where the matrix is
>>>> re-allocated inside SNES iterations so that I can make sure that I do it
>>>> correctly?
>>>>
>>>
>>> No, unfortunately. Contributions always welcome :)
>>>
>>
>>
>> OK, as a test case I'd like to modify snes/tutorials/ex1.c to destroy and
>> reallocate the jacobian matrix with the same sparsity pattern
>> inside FormJacobian1 (once I can do that with the same sparsity pattern,
>> then it should be straightforward to do the same thing with a modified
>> sparsity pattern). To do that, I tried adding the following code inside
>> FormJacobian1:
>>
>>   ierr = MatDestroy(&B);
>>
>
> You do not want to destroy the matrix. There would be no way to get
> another Mat back out of the function.
>
>
>>   ierr = MatCreate(PETSC_COMM_WORLD,&B);CHKERRQ(ierr);
>>
>
> And do not recreate it.
>
>
>>   ierr = MatSetSizes(B,PETSC_DECIDE,PETSC_DECIDE,2,2);CHKERRQ(ierr);
>>
>
> Nor reset the sizes. However, you do want to call
>
>   MatSetPreallocationXAIJ();
>
> which should work
>
>
>>   ierr = MatSetFromOptions(B);CHKERRQ(ierr);
>>   ierr = MatSetUp(B);CHKERRQ(ierr);
>>
>
>   THanks,
>
>     Matt
>

Thanks! Calling MatSetPreallocationXAIJ inside FormJacobian works for me.

So I guess this means that we create a new sparsity pattern "from scratch" each
time we call MatSetPreallocationXAIJ and set the values in the matrix?

Thanks,
David
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20180718/04e5e579/attachment.html>