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

[Matthew Knepley]

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/ccff1d77/attachment-0001.html>