Thanks Matt, for pointing me in the right direction. On Wed, May 19, 2021 at 1:42 PM Matthew Knepley <[email protected]> wrote:
On Wed, May 19, 2021 at 1:30 PM Kaushik Vijaykumar <[email protected]> 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 <[email protected]> wrote:
On Wed, May 19, 2021 at 11:48 AM Kaushik Vijaykumar < [email protected]> 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 <[email protected]> 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 < [email protected]> 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/>