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

[Matthew Knepley]

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/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20210519/766b2a3f/attachment.html>