On Wed, Jul 18, 2018 at 10:15 PM David Knezevic <[email protected]> wrote:
> On Wed, Jul 18, 2018 at 3:34 PM, Matthew Knepley <[email protected]> > wrote: > >> On Wed, Jul 18, 2018 at 3:25 PM David Knezevic < >> [email protected]> wrote: >> >>> On Wed, Jul 18, 2018 at 1:59 PM, Matthew Knepley <[email protected]> >>> wrote: >>> >>>> On Wed, Jul 18, 2018 at 1:31 PM David Knezevic < >>>> [email protected]> 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? > Yes. Matt > Thanks, > David > -- 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.caam.rice.edu/~mk51/>
