[petsc-users] Normalize vectors

[Barry Smith]

On May 11, 2014, at 8:07 PM, LikunTan <tlk0812 at hotmail.com> wrote:

> Dear Petsc developers,
> 
> I have a vector M which consists of a series of 3d vectors, and I want to reset M by normalizing each 3d vector. Here is my code:

   This code is wrong in several ways

> 
> /**********************************************************************
> VecGetArray(M, &aM);     

    This always returns an array that is indexed starting at 0, so the code below when you access with aM[node*3+l] is like totally reading from the wrong place. 

     Instead use 

typedef struct {
   PetscScalar x,y,z;
} Field; 

   Field *aM;

   DMDAVecGetArray(da,M,&aM);

   This routine returns an array whose index starts at xs

> DMDAGetCorners(da, &xs, 0, 0, &xm, 0, 0); 
> 	
> for(node=xs; node<xs+xm; node++)
> {
> 	mag=PetscSqrtScalar(aM[xs].x*aM[xs].x +  aM[xs].y*aM[xs].y  aM[xs].z*aM[xs].z);
>         if (mag != 0.0) {
>            aM[xs].x /= mag;
>            aM[xs].y /= mag;
>            aM[xs].z /= mag;
          }
> }
> DMDAVecRestoreArray(da,M, &aM);       
> 
 When accessing the vector arrays directly you do not need VecAssemblyBegin/End() they are only for use with VecSetValues()

  Also use valgrind to find memory access problems: http://www.mcs.anl.gov/petsc/documentation/faq.html#valgrind

> VecView(M, PETSC_VIEWER_STDOUT_WORLD);
> **********************************************************************/
> 
> but I got the error at the last step:
> --------------------------------------------------------------------------
> mpiexec noticed that process rank 3 with PID 17156 on node compute-21-8.local exited on signal 6 (Aborted).
> --------------------------------------------------------------------------
> and if I commented out VecView, and used vector M for other operations, e.g.
> KSPSolve(ksp, M, b), I got "memory corruption" message. Your comment on this issue is well appreciated.