<div dir="ltr"><div dir="ltr">On Wed, May 19, 2021 at 1:30 PM Kaushik Vijaykumar <<a href="mailto:kaushikv318@gmail.com">kaushikv318@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Matt,<div><br></div><div>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 ?</div></div></blockquote><div><br></div><div>Yes, exactly. You want the residual</div><div><br></div><div> F(U) = 0</div><div><br></div><div>which for you sounds like</div><div><br></div><div> K(U) U - F = 0</div><div><br></div><div>and then the Jacobian callback builds</div><div><br></div><div> K(U)</div><div><br></div><div>from the input U.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Thanks</div><div>Kaushik</div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, May 19, 2021 at 11:51 AM Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div dir="ltr">On Wed, May 19, 2021 at 11:48 AM Kaushik Vijaykumar <<a href="mailto:kaushikv318@gmail.com" target="_blank">kaushikv318@gmail.com</a>> wrote:<br></div><div class="gmail_quote"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr">Barry,<div><br></div><div>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.</div><div>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.</div><div><br></div><div>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().</div><div><br></div><div>Please let me know, if this unclear and I can clarify with more details</div></div></blockquote><div><br></div><div>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.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Thanks</div><div>Kaushik</div><div><br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Wed, May 19, 2021 at 11:25 AM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div><blockquote type="cite"><div dir="ltr"><div><b>ierr = MatSeqAIJGetArray(ptr->K,&Kc);CHKERRQ(ierr); </b></div><div><b>ierr = MatGetLocalSize(ptr->K,&m,&n); CHKERRQ(ierr);<br></b></div><b>for(i=0;i<m;i++)<br> { <br> for(j=0;j<n;j++)<br> {<br> v = *(Kc + (i*n+j));<br> ierr = MatSetValues(Kloc,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);<br> }<br></b><div><b> } </b><br><div> </div></div></div></blockquote><div><br></div>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. <div><br></div><div>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 <a href="https://www.mcs.anl.gov/petsc/documentation/faq.html#functionjacobian" target="_blank">https://www.mcs.anl.gov/petsc/documentation/faq.html#functionjacobian</a> nor do you need to pass into FormFunction in the context the F or Jacobian matrix. </div><div><br></div><div>Maybe if you explained your use case a bit more we could make suggestions on how to accomplish it.</div><div><br></div><div>Barry</div><div><br></div><div><div><br><div><br><blockquote type="cite"><div>On May 19, 2021, at 10:00 AM, Kaushik Vijaykumar <<a href="mailto:kaushikv318@gmail.com" target="_blank">kaushikv318@gmail.com</a>> wrote:</div><br><div><div dir="ltr">Hi everyone,<div><br></div><div>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 "<b>formfunction</b>" to calculate the residual and "<b>formjacobian</b>" to calculate the tangent stiffness matrix. To do so, I have defined a struct:</div><div><br></div><div><b>extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);<br>extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);</b><br><br>// define a struct type to pass K and f as "user context"<br><b>typedef struct<br> {<br> Mat K;<br> Vec f;<br> } K_and_f;</b><br></div><div><br></div><div>In the main program, the struct is declared</div><div><b>int main()</b></div><div><b>{</b></div><div><b>K_and_f main_K_and_f;</b> // declare the struct</div><div><br></div>// SNES - Populate K and f into the struct<br><b>main_K_and_f.K = K;</b> // K matrix<br><div><b>main_K_and_f.f = f; </b>// f vector <br></div><div>....</div><div>....</div><div><b>}</b></div><div>In form function</div><div><br></div><div><b>PetscErrorCode FormFunction(SNES snes, Vec x, Vec F, void* ctx) {<br> PetscErrorCode ierr;<br> PetscReal *ax,*c;<br> PetscReal *aF;<br> PetscScalar *Kc,v; <br> PetscInt nlocal,m,n,i,j,index[100000];<br><br> </b> // Create local Vec and Mat<br><b> Mat Kloc;</b><br><b> Vec Floc, Uloc, KUloc, resloc;</b><br><b> </b><br><b> K_and_f* ptr = (K_and_f*) ctx; </b>// cast the pointer to void into pointer to struct<br></div><div><b> </b>// Get local F array, FLOC</div><b> ierr = VecGhostGetLocalForm(ptr->f,&Floc);CHKERRQ(ierr);</b><div><br></div><div><br></div><div>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. </div><div><br></div><div>The following trial gives me incorrect K values in Kloc:</div><div><b>ierr = MatSeqAIJGetArray(ptr->K,&Kc);CHKERRQ(ierr); </b></div><div><b>ierr = MatGetLocalSize(ptr->K,&m,&n); CHKERRQ(ierr);<br></b></div><b>for(i=0;i<m;i++)<br> { <br> for(j=0;j<n;j++)<br> {<br> v = *(Kc + (i*n+j));<br> ierr = MatSetValues(Kloc,1,&i,1,&j,&v,INSERT_VALUES);CHKERRQ(ierr);<br> }<br></b><div><b> } </b><br><div> </div><div> When I compare K and Kloc, they are not identical </div><div><br></div><div><br></div><div>Please let me know if there is an equivalent function like <b>VecGhostGetLocalForm()</b> for Matrices. If not, is there a better way to do this.</div><div><br></div><div>Any guidance/help is greatly appreciated.</div><div><br></div><div>Thanks</div><div>Kaushik</div><div><br></div><div><br></div><div><br></div><div><br></div></div></div>
</div></blockquote></div><br></div></div></div></blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</blockquote></div>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.cse.buffalo.edu/~knepley/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>