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

Wed May 19 12:50:21 CDT 2021

You may be interested in discussion of computing and representing Newton linearizations for hyperelasticity starting here:

https://libceed.readthedocs.io/en/latest/examples/solids/#id4

Kaushik Vijaykumar <kaushikv318 at gmail.com> writes:

> 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/>
>>