Hello,<div><br></div><div>I am trying to learn SNES to solve some non-linear equations. </div><div><br></div><div>Consider the following non-linear equation which has solution (x,y,z) = (2,3,4)</div><div><br></div><div>x^2 - y^2 + z^2 = 11</div>
<div>xy + y^2 - 3z = 3</div><div>x - xz + yz = 6 </div><div><br></div><div>For solving this, I am able to get the following simple code I wrote to work on 1 processor, but it seems to be failing for 2 or more processors. </div>
<div><br></div><div>The errors I am getting when I run with two processors (mpirun -n 2 ./a.out -snes_monitor )</div><div>Can anybody suggestions to fix this issue? It looks like a very trivial problem, but I am not sure what is wrong.</div>
<div><br></div><div>(The code I wrote is after the error)</div><div><br></div><div>%---------------------------------------------------</div><div>The error is</div><div><div>[0]PETSC ERROR: PetscMallocValidate: error detected at SNESComputeFunction() line 1889 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/interface/snes.c</div>
<div>[0]PETSC ERROR: Memory [id=0(16)] at address 0xa1933f0 is corrupted (probably write past end of array)</div><div>[0]PETSC ERROR: Memory originally allocated in VecCreate_MPI_Private() line 193 in /home/gaurish108/Software/petsc-3.3-p5/src/vec/vec/impls/mpi/pbvec.c</div>
<div>[1]PETSC ERROR: PetscMallocValidate: error detected at SNESComputeFunction() line 1889 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/interface/snes.c</div><div>[1]PETSC ERROR: Memory [id=0(16)] at address 0xa04c350 is corrupted (probably write past end of array)</div>
<div>[1]PETSC ERROR: Memory originally allocated in VecCreate_MPI_Private() line 193 in /home/gaurish108/Software/petsc-3.3-p5/src/vec/vec/impls/mpi/pbvec.c</div><div>[1]PETSC ERROR: --------------------- Error Message ------------------------------------</div>
<div>[1]PETSC ERROR: Memory corruption!</div><div>[1]PETSC ERROR: !</div><div>[1]PETSC ERROR: ------------------------------------------------------------------------</div><div>[1]PETSC ERROR: Petsc Release Version 3.3.0, Patch 5, Sat Dec 1 15:10:41 CST 2012 </div>
<div>[1]PETSC ERROR: See docs/changes/index.html for recent updates.</div><div>[1]PETSC ERROR: See docs/faq.html for hints about trouble shooting.</div><div>[0]PETSC ERROR: --------------------- Error Message ------------------------------------</div>
<div>[0]PETSC ERROR: Memory corruption!</div><div>[0]PETSC ERROR: !</div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: Petsc Release Version 3.3.0, Patch 5, Sat Dec 1 15:10:41 CST 2012 </div>
<div>[0]PETSC ERROR: See docs/changes/index.html for recent updates.</div><div>[0]PETSC ERROR: See docs/faq.html for hints about trouble shooting.</div><div>[0]PETSC ERROR: [1]PETSC ERROR: See docs/index.html for manual pages.</div>
<div>[1]PETSC ERROR: ------------------------------------------------------------------------</div><div>[1]PETSC ERROR: ./gaurish on a linux-gnu named gaurish108 by gaurish108 Fri Mar 22 11:58:01 2013</div><div>[1]PETSC ERROR: Libraries linked from /home/gaurish108/Software/petsc-3.3-p5/linux-gnu-c-yes-debug-refblas/lib</div>
<div>[1]PETSC ERROR: Configure run at Tue Feb 19 20:22:20 2013</div><div>[1]PETSC ERROR: Configure options --download-f-blas-lapack=1</div><div>[1]PETSC ERROR: ------------------------------------------------------------------------</div>
<div>[1]PETSC ERROR: PetscMallocValidate() line 150 in /home/gaurish108/Software/petsc-3.3-p5/src/sys/memory/mtr.c</div><div>[1]PETSC ERROR: SNESComputeFunction() line 1889 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/interface/snes.c</div>
<div>[1]PETSC ERROR: SNESSolve_LS() line 161 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/impls/ls/ls.c</div><div>[1]PETSC ERROR: SNESSolve() line 3536 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/interface/snes.c</div>
<div>[1]PETSC ERROR: main() line 69 in gaurish.c</div><div>See docs/index.html for manual pages.</div><div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: ./gaurish on a linux-gnu named gaurish108 by gaurish108 Fri Mar 22 11:58:01 2013</div>
<div>[0]PETSC ERROR: Libraries linked from /home/gaurish108/Software/petsc-3.3-p5/linux-gnu-c-yes-debug-refblas/lib</div><div>[0]PETSC ERROR: Configure run at Tue Feb 19 20:22:20 2013</div><div>[0]PETSC ERROR: Configure options --download-f-blas-lapack=1</div>
<div>[0]PETSC ERROR: ------------------------------------------------------------------------</div><div>[0]PETSC ERROR: PetscMallocValidate() line 150 in /home/gaurish108/Software/petsc-3.3-p5/src/sys/memory/mtr.c</div><div>
[0]PETSC ERROR: SNESComputeFunction() line 1889 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/interface/snes.c</div><div>[0]PETSC ERROR: SNESSolve_LS() line 161 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/impls/ls/ls.c</div>
<div>[0]PETSC ERROR: SNESSolve() line 3536 in /home/gaurish108/Software/petsc-3.3-p5/src/snes/interface/snes.c</div><div>[0]PETSC ERROR: main() line 69 in gaurish.c</div><div>--------------------------------------------------------------------------</div>
<div>MPI_ABORT was invoked on rank 1 in communicator MPI_COMM_WORLD </div><div>with errorcode 78.</div><div><br></div><div>NOTE: invoking MPI_ABORT causes Open MPI to kill all MPI processes.</div><div>You may or may not see output from other processes, depending on</div>
<div>exactly when Open MPI kills them.</div><div>--------------------------------------------------------------------------</div><div>--------------------------------------------------------------------------</div><div>mpirun has exited due to process rank 1 with PID 17933 on</div>
<div>node gaurish108 exiting without calling "finalize". This may</div><div>have caused other processes in the application to be</div><div>terminated by signals sent by mpirun (as reported here).</div><div>--------------------------------------------------------------------------</div>
<div>[gaurish108:17931] 1 more process has sent help message help-mpi-api.txt / mpi-abort</div><div>[gaurish108:17931] Set MCA parameter "orte_base_help_aggregate" to 0 to see all help / error messages</div></div>
<div><br></div><div><br></div><div><br></div><div>%------------------------------------------------------</div><div>The code is</div><div><br></div><div>static char help[] = "Shows solve a nonlinear equation using Newton iteration method.\n\n";</div>
<div><div>#include <petscsnes.h></div><div>#include <math.h></div><div>#include <petscsys.h></div><div>#include <stdio.h></div><div>#include <petsctime.h></div><div>/* User-defined routines*/</div>
<div>extern PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);</div><div>extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*);</div><div>extern PetscErrorCode FormInitialGuess(Vec);</div><div>extern PetscErrorCode Monitor(SNES,PetscInt,PetscReal,void *);</div>
<div>#undef __FUNCT__</div><div>#define __FUNCT__ "main"</div><div>int main(int argc,char **argv)</div><div>{ SNES snes; /* SNES context */</div><div> Vec x,r,F; /* vectors */</div>
<div> Mat J; /* Jacobian matrix */</div><div> PetscErrorCode ierr;</div><div> PetscInt its,n = 3,i,maxit,maxf;</div><div> PetscMPIInt size;</div><div> PetscScalar iniRHS[3];</div>
<div> </div><div> PetscInitialize(&argc,&argv,(char *)0,help);</div><div> ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr);</div><div> ierr = PetscOptionsGetInt(PETSC_NULL,"-n",&n,PETSC_NULL);CHKERRQ(ierr);</div>
<div> //-------------------------------------------------------------------------------</div><div> ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);</div><div><br></div><div> ierr = VecCreate(PETSC_COMM_WORLD,&x);CHKERRQ(ierr);</div>
<div> ierr = VecSetSizes(x,PETSC_DECIDE,n);CHKERRQ(ierr);</div><div> ierr = VecSetFromOptions(x);CHKERRQ(ierr);</div><div> ierr = VecDuplicate(x,&r);CHKERRQ(ierr);</div><div> ierr = VecDuplicate(x,&F);CHKERRQ(ierr);</div>
<div><br></div><div> ierr = SNESSetFunction(snes,r,FormFunction,(void*)F);CHKERRQ(ierr);</div><div><br></div><div> ierr = MatCreate(PETSC_COMM_WORLD,&J);CHKERRQ(ierr);</div><div> ierr = MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,n,n);CHKERRQ(ierr);</div>
<div> ierr = MatSetFromOptions(J);CHKERRQ(ierr);</div><div> ierr = MatSetUp(J);CHKERRQ(ierr);</div><div> ierr = SNESSetJacobian(snes,J,J,FormJacobian,PETSC_NULL);CHKERRQ(ierr);</div><div><br></div><div> ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);</div>
<div><br></div><div> // Set RHS of equation to F</div><div> iniRHS[0]=11; iniRHS[1]=3; iniRHS[2]=6;</div><div> for (i=0;i<n;i++)</div><div> {</div><div> ierr = VecSetValues(F,1,&i,&iniRHS[i],INSERT_VALUES);CHKERRQ(ierr);</div>
<div> }</div><div> // initial guess</div><div> ierr = FormInitialGuess(x);CHKERRQ(ierr);</div><div><br></div><div> ierr = SNESSolve(snes,PETSC_NULL,x);CHKERRQ(ierr); // The command to call to solve the equation.</div>
<div> </div><div> /* Free work space. All PETSc objects should be destroyed when they are no longer needed*/</div><div> ierr = VecDestroy(&x);CHKERRQ(ierr); ierr = VecDestroy(&r);CHKERRQ(ierr);</div><div> ierr = VecDestroy(&F);CHKERRQ(ierr);</div>
<div> ierr = MatDestroy(&J);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr);</div><div> ierr = PetscFinalize();</div><div> return 0;</div><div>}</div><div><br></div><div><br></div><div>#undef __FUNCT__</div>
<div>#define __FUNCT__ "FormInitialGuess"</div><div>/*</div><div> FormInitialGuess - Computes initial guess.</div><div><br></div><div> Input/Output Parameter:</div><div>. x - the solution vector</div><div>*/</div>
<div>PetscErrorCode FormInitialGuess(Vec x)</div><div>{</div><div> PetscErrorCode ierr;</div><div> PetscScalar pfive =10;</div><div> ierr = VecSet(x,pfive);CHKERRQ(ierr);</div><div> return 0;</div><div>}</div><div>
<br></div><div><br></div><div>/* ------------------------------------------------------------------- */</div><div>#undef __FUNCT__</div><div>#define __FUNCT__ "FormFunction"</div><div>/* </div><div> FormFunction - Evaluates nonlinear function, F(x).</div>
<div><br></div><div> Input Parameters:</div><div>. snes - the SNES context</div><div>. x - input vector</div><div>. ctx - optional user-defined context, as set by SNESSetFunction()</div><div><br></div><div> Output Parameter:</div>
<div>. f - function vector</div><div><br></div><div> Note:</div><div> The user-defined context can contain any application-specific data</div><div> needed for the function evaluation (such as various parameters, work</div>
<div> vectors, and grid information). In this program the context is just</div><div> a vector containing the right-hand-side of the discretized PDE.</div><div> */</div><div><br></div><div>PetscErrorCode FormFunction(SNES snes,Vec x,Vec f,void *ctx)</div>
<div>{</div><div> Vec g = (Vec)ctx;</div><div> PetscScalar *xx,*ff,*gg;</div><div> PetscErrorCode ierr;</div><div> PetscInt n;</div><div><br></div><div> /*</div><div> Get pointers to vector data.</div>
<div> - For default PETSc vectors, VecGetArray() returns a pointer to</div><div> the data array. Otherwise, the routine is implementation dependent.</div><div> - You MUST call VecRestoreArray() when you no longer need access to</div>
<div> the array.</div><div> */</div><div> ierr = VecGetArray(x,&xx);CHKERRQ(ierr);</div><div> ierr = VecGetArray(f,&ff);CHKERRQ(ierr);</div><div> ierr = VecGetArray(g,&gg);CHKERRQ(ierr);</div><div>
<br></div><div> /*</div><div> Compute function</div><div> */</div><div> ierr = VecGetSize(x,&n);CHKERRQ(ierr);</div><div><br></div><div> ff[0]=pow(xx[0],2)-pow(xx[1],2)+pow(xx[2],2)-gg[0];</div><div> ff[1]=xx[0]*xx[1]+pow(xx[1],2)-3.0*xx[2]-gg[1];</div>
<div> ff[2]=xx[0]-xx[0]*xx[2]+xx[1]*xx[2]-gg[2];</div><div><br></div><div> /*</div><div> Restore vectors</div><div> */</div><div> ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);</div><div> ierr = VecRestoreArray(f,&ff);CHKERRQ(ierr);</div>
<div> ierr = VecRestoreArray(g,&gg);CHKERRQ(ierr);</div><div> return 0;</div><div>}</div><div><br></div><div><br></div><div><br></div><div>/* ------------------------------------------------------------------- */</div>
<div>#undef __FUNCT__</div><div>#define __FUNCT__ "FormJacobian"</div><div>/*</div><div> FormJacobian - Evaluates Jacobian matrix.</div><div><br></div><div> Input Parameters:</div><div>. snes - the SNES context</div>
<div>. x - input vector</div><div>. dummy - optional user-defined context (not used here)</div><div><br></div><div> Output Parameters:</div><div>. jac - Jacobian matrix</div><div>. B - optionally different preconditioning matrix</div>
<div>. flag - flag indicating matrix structure</div><div>*/</div><div><br></div><div>PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *jac,Mat *B,MatStructure*flag,void *dummy)</div><div>{</div><div> PetscScalar *xx,A[3];</div>
<div> PetscErrorCode ierr;</div><div> PetscInt i,n,j[3];</div><div><br></div><div> /*</div><div> Get pointer to vector data</div><div> */</div><div> ierr = VecGetArray(x,&xx);CHKERRQ(ierr);</div><div><br>
</div><div> /*</div><div> Compute Jacobian entries and insert into matrix.</div><div> - Note that in this case we set all elements for a particular</div><div> row at once.</div><div> */</div><div> ierr = VecGetSize(x,&n);CHKERRQ(ierr);</div>
<div> j[0] = 0; j[1] = 1; j[2] = 2; </div><div> /*set the first row of jacobian*/</div><div> i=0;</div><div> A[0]=2*xx[0]; A[1]=-2*xx[1]; A[2]=2*xx[2];</div><div> ierr = MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr);</div>
<div> /*set the second row of jacobian*/</div><div> i=1;</div><div> A[0]=xx[1]; A[1]=xx[0]+2*xx[1]; A[2]=-3.0;</div><div> ierr = MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr);</div><div> /*set the third row of jacobian*/</div>
<div> i=2;</div><div> A[0]=1-xx[2]; A[1]=xx[2]; A[2]=-xx[0]+xx[1];</div><div> ierr = MatSetValues(*jac,1,&i,3,j,A,INSERT_VALUES);CHKERRQ(ierr);</div><div> </div><div><br></div><div><br></div><div> ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr);</div>
<div><br></div><div> /*</div><div> Assemble matrix</div><div> */</div><div> ierr = MatAssemblyBegin(*jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</div><div> ierr = MatAssemblyEnd(*jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr);</div>
<div><br></div><div> return 0;</div><div>}</div><div><br></div><div><br></div></div><div><br></div><div><br></div><div><br></div><div><br></div>