[petsc-users] Matrix assembly in SNES problem

[Barry Smith]

> On Jul 22, 2020, at 7:03 AM, Adolfo Rodriguez <adantra at gmail.com> wrote:
> 
> I am trying to replace the non-linear solver in a flow simulation problem where the matrix sparsity can change during the iterations. I tried successfully to create the matrix within the FormJacobian function but I have a memory leak which I don't know how to fix.
> 
> My FormJacobian function is something like this:
> 
> .. FormJacobian( ..){
>      matrixconstruction(); // returns csr matrix in c++ regular vectors

      What do you mean by c++ regular vectors? Do you mean basic arrays you obtain with New or malloc() or do you mean some std: vector objects? I am guessing just arrays obtained with new.

>      MatSeqAIJSetPreallocationCSR ...;
>      //loop over the entries and 
>      MatSetValues;
>      MatAssemblyBegin();
>      MatAssemblyEnd();
>      return 0;
> }

1) side note

    The routine MatSeqAIJSetPreallocationCSR() which takes ii, jj,vv arguments not only preallocates the memory but also fills up the matrix so you don't need the part of the code
>     //loop over the entries and 
>      MatSetValues;
>      MatAssemblyBegin();
>      MatAssemblyEnd();
    (Yes the function name is confusing)
    
    You should remove it and still get the same results but faster.

2) the memory leak could be in your code or in PETSc. Are you freeing the arrays you get from matrixconstruction? 

  You can run a smallish problem with -malloc_debug and at PetscFinalize() it will print all the memory PETSc has allocated and not freed, and where in the code it was allocated. If the leak is in PETSc or from not freeing PETSc objects then this should find the problem. If the leak is in your C++ allocations you will need to use the Microsoft tools.

  Good luck,


   Barry





> 
> This approach works for small problems but, when dealing with large ones I observe a memory leak. I am pretty sure that It is related to the fact that I am creating new matrices every time FormJacobian is called but I don't know how to destroy it between iterations.
> 
> Any suggestions?
> 
> Regards,
> 
> Adolfo