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 ?
 

Thanks
Kaushik

On Wed, May 19, 2021 at 11:51 AM Matthew Knepley <knepley@gmail.com> wrote:
On Wed, May 19, 2021 at 11:48 AM Kaushik Vijaykumar <kaushikv318@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@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@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