[petsc-users] Using SNES to solve Non linear elasticity, passing matrix to form function

Kaushik Vijaykumar kaushikv318 at gmail.com
Wed May 19 12:46:15 CDT 2021 

Thanks Matt, for pointing me in the right direction.



On Wed, May 19, 2021 at 1:42 PM Matthew Knepley <knepley at gmail.com> wrote:

> On Wed, May 19, 2021 at 1:30 PM Kaushik Vijaykumar <kaushikv318 at gmail.com>
> wrote:
>
>> Matt,
>>
>> Thanks for the comment, I see your point now. I was using the linear
>> elastic [K] from the previous load step as the guess for the jacobian. The
>> stiffness matrix was being built outside formfunction using displacements
>> from the previous load step. What is the best approach, Build [K] abd {F}
>> inside the form function and compute the residual [K]_{i}{U}_{i} - {F}_{i},
>> where i designates the iterated matrices and vector for a given load step ?
>>
>
> Yes, exactly. You want the residual
>
>   F(U) = 0
>
> which for you sounds like
>
>   K(U) U - F = 0
>
> and then the Jacobian callback builds
>
>   K(U)
>
> from the input U.
>
>   Thanks,
>
>      Matt
>
>
>> Thanks
>> Kaushik
>>
>> On Wed, May 19, 2021 at 11:51 AM Matthew Knepley <knepley at gmail.com>
>> wrote:
>>
>>> On Wed, May 19, 2021 at 11:48 AM Kaushik Vijaykumar <
>>> kaushikv318 at gmail.com> wrote:
>>>
>>>> Barry,
>>>>
>>>> Thanks for your response. I am trying to copy K into Kloc, the purpose
>>>> of doing that is to calculate the residual in formfunction - Residual =
>>>> [Kloc]{X} - {Floc}, where Kloc and Floc are the copies of assembled global
>>>> K matrix and F vector, and X is the current iterated solution. This is the
>>>> reason that I need to access the global stiffness matrix and force vector
>>>> in formfunction.
>>>> I agree that formjacobian only needs access to global stiffness matrix
>>>> K to form the jacobian and does not need to access the force vector.
>>>>
>>>> Therefore, to be able to pass both K and F to formfunction, I defined a
>>>> struct that contains a vector F and matrix K and populated them in the
>>>> main().
>>>>
>>>> Please let me know, if this unclear and I can clarify with more details
>>>>
>>>
>>> This still does not make sense. If you have a nonlinear problem, how can
>>> you know the Jacobian up front? It changes with the input point.
>>>
>>>   Thanks,
>>>
>>>     Matt
>>>
>>>
>>>> Thanks
>>>> Kaushik
>>>>
>>>>
>>>> On Wed, May 19, 2021 at 11:25 AM Barry Smith <bsmith at petsc.dev> wrote:
>>>>
>>>>> *ierr = MatSeqAIJGetArray(ptr->K,&Kc);CHKERRQ(ierr); *
>>>>>
>>>>> *ierr = MatGetLocalSize(ptr->K,&m,&n); CHKERRQ(ierr);*
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> *for(i=0;i<m;i++)    {       for(j=0;j<n;j++)      {        v  = *(Kc
>>>>> + (i*n+j));        ierr =
>>>>> MatSetValues(Kloc,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);      }*
>>>>> *    } *
>>>>>
>>>>>
>>>>>
>>>>> I don't understand the purpose of this code fragment. Are you wanting
>>>>> to copy your K into Kloc? What is the purpose of Kloc? With standard usage
>>>>> of SNES one only needs to provide PETSc a global matrix, one does not need
>>>>> to work with "ghosted" matrices at all.
>>>>>
>>>>> Also normally one fills the f vector in FormFunction and the Jacobian
>>>>> in FormJacobian and there is no need to access the Jacobians at all in SNES
>>>>> https://www.mcs.anl.gov/petsc/documentation/faq.html#functionjacobian
>>>>> nor do you need to pass into FormFunction in the context the F or Jacobian
>>>>> matrix.
>>>>>
>>>>> Maybe if you explained your use case a bit more we could make
>>>>> suggestions on how to accomplish it.
>>>>>
>>>>> Barry
>>>>>
>>>>>
>>>>>
>>>>> On May 19, 2021, at 10:00 AM, Kaushik Vijaykumar <
>>>>> kaushikv318 at gmail.com> wrote:
>>>>>
>>>>> Hi everyone,
>>>>>
>>>>> We are in the process of integrating nonlinear solution to our FE code
>>>>> through SNES. The first steps that I understood that need to be done is to
>>>>> be able to pass the assembled stiffness matrix K and force vector F to the "
>>>>> *formfunction*" to calculate the residual and "*formjacobian*" to
>>>>> calculate the tangent stiffness matrix. To do so, I have defined a struct:
>>>>>
>>>>>
>>>>> *extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);extern
>>>>> PetscErrorCode FormFunction(SNES,Vec,Vec,void*);*
>>>>>
>>>>> // define a struct type to pass K and f as "user context"
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> *typedef struct    {      Mat K;      Vec f;    } K_and_f;*
>>>>>
>>>>> In the main program, the struct is declared
>>>>> *int main()*
>>>>> *{*
>>>>> *K_and_f main_K_and_f;* // declare the struct
>>>>>
>>>>> // SNES - Populate K and f into the struct
>>>>> *main_K_and_f.K = K;* // K matrix
>>>>> *main_K_and_f.f = f; *// f vector
>>>>> ....
>>>>> ....
>>>>> *}*
>>>>> In form function
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> *PetscErrorCode FormFunction(SNES snes, Vec x, Vec F, void* ctx) {
>>>>> PetscErrorCode ierr;    PetscReal *ax,*c;    PetscReal *aF;    PetscScalar
>>>>> *Kc,v;     PetscInt nlocal,m,n,i,j,index[100000];   * // Create local
>>>>> Vec and Mat
>>>>> *    Mat Kloc;*
>>>>> *    Vec Floc, Uloc, KUloc, resloc;*
>>>>>
>>>>> *    K_and_f* ptr = (K_and_f*) ctx; *// cast the pointer to void into
>>>>> pointer to struct
>>>>>     // Get local F array, FLOC
>>>>> *    ierr = VecGhostGetLocalForm(ptr->f,&Floc);CHKERRQ(ierr);*
>>>>>
>>>>>
>>>>> I am able to get the f array from the main program using "
>>>>> VecGhostGetLocalForm" and the vec is correct. However, I am having trouble
>>>>> with reading the K matrix in formfunction.
>>>>>
>>>>> The following trial gives me incorrect K values in Kloc:
>>>>> *ierr = MatSeqAIJGetArray(ptr->K,&Kc);CHKERRQ(ierr); *
>>>>>
>>>>> *ierr = MatGetLocalSize(ptr->K,&m,&n); CHKERRQ(ierr);*
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>>
>>>>> *for(i=0;i<m;i++)    {       for(j=0;j<n;j++)      {        v  = *(Kc
>>>>> + (i*n+j));        ierr =
>>>>> MatSetValues(Kloc,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);      }*
>>>>> *    } *
>>>>>
>>>>>     When I compare K and Kloc, they are not identical
>>>>>
>>>>>
>>>>> Please let me know if there is an equivalent function like
>>>>> *VecGhostGetLocalForm()* for Matrices. If not, is there a better way
>>>>> to do this.
>>>>>
>>>>> Any guidance/help is greatly appreciated.
>>>>>
>>>>> 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/>
>>>
>>
>
> --
> 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/20210519/f89445cc/attachment.html>

More information about the petsc-users mailing list