Actual source code: ex10.c

  1: /* 
  2:   Program usage:  mpiexec -n <procs> usg [-help] [all PETSc options] 
  3: */

  5: #if !defined(PETSC_USE_COMPLEX)

  7: static char help[] = "An Unstructured Grid Example.\n\
  8: This example demonstrates how to solve a nonlinear system in parallel\n\
  9: with SNES for an unstructured mesh. The mesh and partitioning information\n\
 10: is read in an application defined ordering,which is later transformed\n\
 11: into another convenient ordering (called the local ordering). The local\n\
 12: ordering, apart from being efficient on cpu cycles and memory, allows\n\
 13: the use of the SPMD model of parallel programming. After partitioning\n\
 14: is done, scatters are created between local (sequential)and global\n\
 15: (distributed) vectors. Finally, we set up the nonlinear solver context\n\
 16: in the usual way as a structured grid  (see\n\
 17: petsc/src/snes/examples/tutorials/ex5.c).\n\
 18: This example also illustrates the use of parallel matrix coloring.\n\
 19: The command line options include:\n\
 20:   -vert <Nv>, where Nv is the global number of nodes\n\
 21:   -elem <Ne>, where Ne is the global number of elements\n\
 22:   -nl_par <lambda>, where lambda is the multiplier for the non linear term (u*u) term\n\
 23:   -lin_par <alpha>, where alpha is the multiplier for the linear term (u)\n\
 24:   -fd_jacobian_coloring -mat_coloring_type lf\n";

 26: /*T
 27:    Concepts: SNES^unstructured grid
 28:    Concepts: AO^application to PETSc ordering or vice versa;
 29:    Concepts: VecScatter^using vector scatter operations;
 30:    Processors: n
 31: T*/

 33: /* ------------------------------------------------------------------------

 35:    PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian.

 37:    The Laplacian is approximated in the following way: each edge is given a weight
 38:    of one meaning that the diagonal term will have the weight equal to the degree
 39:    of a node. The off diagonal terms will get a weight of -1. 

 41:    -----------------------------------------------------------------------*/

 43: /*
 44:    Include petscao.h so that we can use AO (Application Ordering) object's services.
 45:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
 46:    file automatically includes:
 47:      petscsys.h       - base PETSc routines   petscvec.h - vectors
 48:      petscmat.h - matrices
 49:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
 50:      petscviewer.h - viewers               petscpc.h  - preconditioners
 51:      petscksp.h   - linear solvers
 52: */
 53: #include "petscao.h"
 54: #include "petscsnes.h"


 57: #define MAX_ELEM      500  /* Maximum number of elements */
 58: #define MAX_VERT      100  /* Maximum number of vertices */
 59: #define MAX_VERT_ELEM   3  /* Vertices per element       */

 61: /*
 62:   Application-defined context for problem specific data 
 63: */
 64: typedef struct {
 65:       PetscInt    Nvglobal,Nvlocal;            /* global and local number of vertices */
 66:       PetscInt    Neglobal,Nelocal;            /* global and local number of vertices */
 67:       PetscInt           AdjM[MAX_VERT][50];           /* adjacency list of a vertex */
 68:       PetscInt           itot[MAX_VERT];               /* total number of neighbors for a vertex */
 69:       PetscInt           icv[MAX_ELEM][MAX_VERT_ELEM]; /* vertices belonging to an element */
 70:       PetscInt          v2p[MAX_VERT];                /* processor number for a vertex */
 71:       PetscInt    *locInd,*gloInd;             /* local and global orderings for a node */
 72:       Vec           localX,localF;               /* local solution (u) and f(u) vectors */
 73:       PetscReal          non_lin_param;                /* nonlinear parameter for the PDE */
 74:       PetscReal          lin_param;                    /* linear parameter for the PDE */
 75:       VecScatter  scatter;                      /* scatter context for the local and 
 76:                                                     distributed vectors */
 77: } AppCtx;

 79: /*
 80:   User-defined routines
 81: */
 82: PetscErrorCode  FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
 83:                 FormFunction(SNES,Vec,Vec,void*),
 84:                 FormInitialGuess(AppCtx*,Vec);

 88: int main(int argc,char **argv)
 89: {
 90:   SNES                   snes;                 /* SNES context */
 91:   const SNESType         type = SNESLS;        /* default nonlinear solution method */
 92:   Vec                    x,r;                  /* solution, residual vectors */
 93:   Mat                    Jac;                  /* Jacobian matrix */
 94:   AppCtx                 user;                 /* user-defined application context */
 95:   AO                     ao;                   /* Application Ordering object */
 96:   IS                     isglobal,islocal;     /* global and local index sets */
 97:   PetscMPIInt            rank,size;            /* rank of a process, number of processors */
 98:   PetscInt               rstart;               /* starting index of PETSc ordering for a processor */
 99:   PetscInt               nfails;               /* number of unsuccessful Newton steps */
100:   PetscInt               bs = 1;               /* block size for multicomponent systems */
101:   PetscInt               nvertices;            /* number of local plus ghost nodes of a processor */
102:   PetscInt               *pordering;           /* PETSc ordering */
103:   PetscInt               *vertices;            /* list of all vertices (incl. ghost ones) 
104:                                                 on a processor */
105:   PetscInt               *verticesmask,*svertices;
106:   PetscInt               *tmp;
107:   PetscInt               i,j,jstart,inode,nb,nbrs,Nvneighborstotal = 0;
108:   PetscErrorCode         ierr;
109:   PetscInt               its,N;
110:   PetscScalar            *xx;
111:   char                   str[256],form[256],part_name[256];
112:   FILE                   *fptr,*fptr1;
113:   ISLocalToGlobalMapping isl2g;
114:   int                    dtmp;
115: #if defined (UNUSED_VARIABLES)
116:   PetscDraw              draw;                 /* drawing context */
117:   PetscScalar            *ff,*gg;
118:   PetscReal              tiny = 1.0e-10,zero = 0.0,one = 1.0,big = 1.0e+10;
119:   PetscInt               *tmp1,*tmp2;
120: #endif
121:   MatFDColoring          matfdcoloring = 0;
122:   PetscTruth             fd_jacobian_coloring = PETSC_FALSE;

124:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
125:      Initialize program
126:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

128:   PetscInitialize(&argc,&argv,"options.inf",help);
129:   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
130:   MPI_Comm_size(MPI_COMM_WORLD,&size);

132:   /* The current input file options.inf is for 2 proc run only */
133:   if (size != 2) SETERRQ(1,"This Example currently runs on 2 procs only.");

135:   /*
136:      Initialize problem parameters
137:   */
138:   user.Nvglobal = 16;      /*Global # of vertices  */
139:   user.Neglobal = 18;      /*Global # of elements  */
140:   PetscOptionsGetInt(PETSC_NULL,"-vert",&user.Nvglobal,PETSC_NULL);
141:   PetscOptionsGetInt(PETSC_NULL,"-elem",&user.Neglobal,PETSC_NULL);
142:   user.non_lin_param = 0.06;
143:   PetscOptionsGetReal(PETSC_NULL,"-nl_par",&user.non_lin_param,PETSC_NULL);
144:   user.lin_param = -1.0;
145:   PetscOptionsGetReal(PETSC_NULL,"-lin_par",&user.lin_param,PETSC_NULL);
146:   user.Nvlocal = 0;
147:   user.Nelocal = 0;

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:       Read the mesh and partitioning information  
151:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: 
153:   /*
154:      Read the mesh and partitioning information from 'adj.in'.
155:      The file format is as follows.
156:      For each line the first entry is the processor rank where the 
157:      current node belongs. The second entry is the number of 
158:      neighbors of a node. The rest of the line is the adjacency
159:      list of a node. Currently this file is set up to work on two 
160:      processors.

162:      This is not a very good example of reading input. In the future, 
163:      we will put an example that shows the style that should be 
164:      used in a real application, where partitioning will be done 
165:      dynamically by calling partitioning routines (at present, we have 
166:      a  ready interface to ParMeTiS). 
167:    */
168:   fptr = fopen("adj.in","r");
169:   if (!fptr) {
170:       SETERRQ(0,"Could not open adj.in")
171:   }
172: 
173:   /*
174:      Each processor writes to the file output.<rank> where rank is the
175:      processor's rank.
176:   */
177:   sprintf(part_name,"output.%d",rank);
178:   fptr1 = fopen(part_name,"w");
179:   if (!fptr1) {
180:       SETERRQ(0,"Could no open output file");
181:   }
182:   PetscMalloc(user.Nvglobal*sizeof(PetscInt),&user.gloInd);
183:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Rank is %D\n",rank);
184:   for (inode = 0; inode < user.Nvglobal; inode++) {
185:     fgets(str,256,fptr);
186:     sscanf(str,"%d",&dtmp);user.v2p[inode] = dtmp;
187:     if (user.v2p[inode] == rank) {
188:        PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %D belongs to processor %D\n",inode,user.v2p[inode]);
189:        user.gloInd[user.Nvlocal] = inode;
190:        sscanf(str,"%*d %d",&dtmp);nbrs = dtmp;
191:        PetscFPrintf(PETSC_COMM_SELF,fptr1,"Number of neighbors for the vertex %D is %D\n",inode,nbrs);
192:        user.itot[user.Nvlocal] = nbrs;
193:        Nvneighborstotal += nbrs;
194:        for (i = 0; i < user.itot[user.Nvlocal]; i++){
195:          form[0]='\0';
196:          for (j=0; j < i+2; j++){
197:            PetscStrcat(form,"%*d ");
198:          }
199:            PetscStrcat(form,"%d");
200:            sscanf(str,form,&dtmp); user.AdjM[user.Nvlocal][i] = dtmp;
201:            PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[user.Nvlocal][i]);
202:         }
203:         PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
204:         user.Nvlocal++;
205:      }
206:    }
207:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Total # of Local Vertices is %D \n",user.Nvlocal);

209:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
210:      Create different orderings
211:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

213:   /*
214:     Create the local ordering list for vertices. First a list using the PETSc global 
215:     ordering is created. Then we use the AO object to get the PETSc-to-application and 
216:     application-to-PETSc mappings. Each vertex also gets a local index (stored in the 
217:     locInd array). 
218:   */
219:   MPI_Scan(&user.Nvlocal,&rstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
220:   rstart -= user.Nvlocal;
221:   PetscMalloc(user.Nvlocal*sizeof(PetscInt),&pordering);

223:   for (i=0; i < user.Nvlocal; i++) {
224:     pordering[i] = rstart + i;
225:   }

227:   /* 
228:     Create the AO object 
229:   */
230:   AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
231:   PetscFree(pordering);
232: 
233:   /* 
234:     Keep the global indices for later use 
235:   */
236:   PetscMalloc(user.Nvlocal*sizeof(PetscInt),&user.locInd);
237:   PetscMalloc(Nvneighborstotal*sizeof(PetscInt),&tmp);
238: 
239:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
240:     Demonstrate the use of AO functionality 
241:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

243:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Before AOApplicationToPetsc, local indices are : \n");
244:   for (i=0; i < user.Nvlocal; i++) {
245:    PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.gloInd[i]);
246:    user.locInd[i] = user.gloInd[i];
247:   }
248:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
249:   jstart = 0;
250:   for (i=0; i < user.Nvlocal; i++) {
251:    PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.gloInd[i]);
252:    for (j=0; j < user.itot[i]; j++) {
253:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
254:     tmp[j + jstart] = user.AdjM[i][j];
255:    }
256:    jstart += user.itot[i];
257:    PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
258:   }

260:   /* 
261:     Now map the vlocal and neighbor lists to the PETSc ordering 
262:   */
263:   AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);
264:   AOApplicationToPetsc(ao,Nvneighborstotal,tmp);
265:   AODestroy(ao);

267:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"After AOApplicationToPetsc, local indices are : \n");
268:   for (i=0; i < user.Nvlocal; i++) {
269:    PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.locInd[i]);
270:   }
271:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");

273:   jstart = 0;
274:   for (i=0; i < user.Nvlocal; i++) {
275:    PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.locInd[i]);
276:    for (j=0; j < user.itot[i]; j++) {
277:     user.AdjM[i][j] = tmp[j+jstart];
278:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
279:    }
280:    jstart += user.itot[i];
281:    PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
282:   }

284:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
285:      Extract the ghost vertex information for each processor
286:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
287:   /*
288:    Next, we need to generate a list of vertices required for this processor
289:    and a local numbering scheme for all vertices required on this processor.
290:       vertices - integer array of all vertices needed on this processor in PETSc 
291:                  global numbering; this list consists of first the "locally owned" 
292:                  vertices followed by the ghost vertices.
293:       verticesmask - integer array that for each global vertex lists its local 
294:                      vertex number (in vertices) + 1. If the global vertex is not
295:                      represented on this processor, then the corresponding
296:                      entry in verticesmask is zero
297:  
298:       Note: vertices and verticesmask are both Nvglobal in length; this may
299:     sound terribly non-scalable, but in fact is not so bad for a reasonable
300:     number of processors. Importantly, it allows us to use NO SEARCHING
301:     in setting up the data structures.
302:   */
303:   PetscMalloc(user.Nvglobal*sizeof(PetscInt),&vertices);
304:   PetscMalloc(user.Nvglobal*sizeof(PetscInt),&verticesmask);
305:   PetscMemzero(verticesmask,user.Nvglobal*sizeof(PetscInt));
306:   nvertices = 0;
307: 
308:   /* 
309:     First load "owned vertices" into list 
310:   */
311:   for (i=0; i < user.Nvlocal; i++) {
312:     vertices[nvertices++]   = user.locInd[i];
313:     verticesmask[user.locInd[i]] = nvertices;
314:   }
315: 
316:   /* 
317:     Now load ghost vertices into list 
318:   */
319:   for (i=0; i < user.Nvlocal; i++) {
320:     for (j=0; j < user.itot[i]; j++) {
321:       nb = user.AdjM[i][j];
322:       if (!verticesmask[nb]) {
323:         vertices[nvertices++] = nb;
324:         verticesmask[nb]      = nvertices;
325:       }
326:     }
327:   }

329:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
330:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n");
331:   for (i=0; i < nvertices; i++) {
332:    PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",vertices[i]);
333:    }
334:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
335: 
336:   /*
337:      Map the vertices listed in the neighbors to the local numbering from
338:     the global ordering that they contained initially.
339:   */
340:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
341:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"After mapping neighbors in the local contiguous ordering\n");
342:   for (i=0; i<user.Nvlocal; i++) {
343:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are :\n",i);
344:     for (j = 0; j < user.itot[i]; j++) {
345:       nb = user.AdjM[i][j];
346:       user.AdjM[i][j] = verticesmask[nb] - 1;
347:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
348:     }
349:    PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
350:   }

352:   N = user.Nvglobal;
353: 
354:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
355:      Create vector and matrix data structures
356:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

358:   /* 
359:     Create vector data structures 
360:   */
361:   VecCreate(MPI_COMM_WORLD,&x);
362:   VecSetSizes(x,user.Nvlocal,N);
363:   VecSetFromOptions(x);
364:   VecDuplicate(x,&r);
365:   VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
366:   VecDuplicate(user.localX,&user.localF);

368:   /*
369:     Create the scatter between the global representation and the 
370:     local representation
371:   */
372:   ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
373:   PetscMalloc(nvertices*sizeof(PetscInt),&svertices);
374:   for (i=0; i<nvertices; i++) svertices[i] = bs*vertices[i];
375:   ISCreateBlock(MPI_COMM_SELF,bs,nvertices,svertices,&isglobal);
376:   PetscFree(svertices);
377:   VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);
378:   ISDestroy(isglobal);
379:   ISDestroy(islocal);

381:   /* 
382:      Create matrix data structure; Just to keep the example simple, we have not done any 
383:      preallocation of memory for the matrix. In real application code with big matrices,
384:      preallocation should always be done to expedite the matrix creation. 
385:   */
386:   MatCreate(MPI_COMM_WORLD,&Jac);
387:   MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
388:   MatSetFromOptions(Jac);

390:   /* 
391:     The following routine allows us to set the matrix values in local ordering 
392:   */
393:   ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs*nvertices,vertices,&isl2g);
394:   MatSetLocalToGlobalMapping(Jac,isl2g);

396:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
397:      Create nonlinear solver context
398:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

400:   SNESCreate(MPI_COMM_WORLD,&snes);
401:   SNESSetType(snes,type);

403:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
404:     Set routines for function and Jacobian evaluation
405:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
406:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

408:    PetscOptionsGetTruth(PETSC_NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,0);
409:    if (!fd_jacobian_coloring){
410:      SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user);
411:    } else { /* Use matfdcoloring */
412:      ISColoring    iscoloring;
413:      MatStructure  flag;
414:      /* Get the data structure of Jac */
415:      FormJacobian(snes,x,&Jac,&Jac,&flag,&user);
416:      /* Create coloring context */
417:      MatGetColoring(Jac,MATCOLORING_SL,&iscoloring);
418:      MatFDColoringCreate(Jac,iscoloring,&matfdcoloring);
419:      MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
420:      MatFDColoringSetFromOptions(matfdcoloring);
421:      /* MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD); */
422:      SNESSetJacobian(snes,Jac,Jac,SNESDefaultComputeJacobianColor,matfdcoloring);
423:      ISColoringDestroy(iscoloring);
424:    }

426:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
427:      Customize nonlinear solver; set runtime options
428:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

430:   SNESSetFromOptions(snes);

432:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
433:      Evaluate initial guess; then solve nonlinear system
434:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

436:   /*
437:      Note: The user should initialize the vector, x, with the initial guess
438:      for the nonlinear solver prior to calling SNESSolve().  In particular,
439:      to employ an initial guess of zero, the user should explicitly set
440:      this vector to zero by calling VecSet().
441:   */
442:   FormInitialGuess(&user,x);

444:    /* 
445:      Print the initial guess 
446:    */
447:   VecGetArray(x,&xx);
448:   for (inode = 0; inode < user.Nvlocal; inode++){
449:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Initial Solution at node %D is %f \n",inode,xx[inode]);
450:   }
451:   VecRestoreArray(x,&xx);

453:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
454:      Now solve the nonlinear system
455:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

457:   SNESSolve(snes,PETSC_NULL,x);
458:   SNESGetIterationNumber(snes,&its);
459:   SNESGetNonlinearStepFailures(snes,&nfails);
460: 
461:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
462:     Print the output : solution vector and other information
463:      Each processor writes to the file output.<rank> where rank is the
464:      processor's rank.
465:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

467:   VecGetArray(x,&xx);
468:   for (inode = 0; inode < user.Nvlocal; inode++)
469:    PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %D is %f \n",inode,xx[inode]);
470:   VecRestoreArray(x,&xx);
471:   fclose(fptr1);
472:   PetscPrintf(MPI_COMM_WORLD,"number of Newton iterations = %D, ",its);
473:   PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %D\n",nfails);

475:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
476:      Free work space.  All PETSc objects should be destroyed when they
477:      are no longer needed.
478:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
479:   PetscFree(user.gloInd);
480:   PetscFree(user.locInd);
481:   PetscFree(vertices);
482:   PetscFree(verticesmask);
483:   PetscFree(tmp);
484:   VecScatterDestroy(user.scatter);
485:   ISLocalToGlobalMappingDestroy(isl2g);
486:   VecDestroy(x);
487:   VecDestroy(r);
488:   VecDestroy(user.localX);
489:   VecDestroy(user.localF);
490:   MatDestroy(Jac);  SNESDestroy(snes);
491:   /*PetscDrawDestroy(draw);*/
492:   if (fd_jacobian_coloring){
493:     MatFDColoringDestroy(matfdcoloring);
494:   }
495:   PetscFinalize();

497:   return 0;
498: }
501: /* --------------------  Form initial approximation ----------------- */

503: /* 
504:    FormInitialGuess - Forms initial approximation.

506:    Input Parameters:
507:    user - user-defined application context
508:    X - vector

510:    Output Parameter:
511:    X - vector
512:  */
513: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
514: {
515:   PetscInt    i,Nvlocal,ierr;
516:   PetscInt    *gloInd;
517:   PetscScalar *x;
518: #if defined (UNUSED_VARIABLES)
519:   PetscReal  temp1,temp,hx,hy,hxdhy,hydhx,sc;
520:   PetscInt   Neglobal,Nvglobal,j,row;
521:   PetscReal  alpha,lambda;

523:   Nvglobal = user->Nvglobal;
524:   Neglobal = user->Neglobal;
525:   lambda   = user->non_lin_param;
526:   alpha    = user->lin_param;
527: #endif

529:   Nvlocal  = user->Nvlocal;
530:   gloInd   = user->gloInd;

532:   /*
533:      Get a pointer to vector data.
534:        - For default PETSc vectors, VecGetArray() returns a pointer to
535:          the data array.  Otherwise, the routine is implementation dependent.
536:        - You MUST call VecRestoreArray() when you no longer need access to
537:          the array.
538:   */
539:   VecGetArray(X,&x);

541:   /*
542:      Compute initial guess over the locally owned part of the grid
543:   */
544:   for (i=0; i < Nvlocal; i++) {
545:     x[i] = (PetscReal)gloInd[i];
546:   }

548:   /*
549:      Restore vector
550:   */
551:   VecRestoreArray(X,&x);
552:   return 0;
553: }
556: /* --------------------  Evaluate Function F(x) --------------------- */
557: /* 
558:    FormFunction - Evaluates nonlinear function, F(x).

560:    Input Parameters:
561: .  snes - the SNES context
562: .  X - input vector
563: .  ptr - optional user-defined context, as set by SNESSetFunction()

565:    Output Parameter:
566: .  F - function vector
567:  */
568: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
569: {
570:   AppCtx         *user = (AppCtx*)ptr;
572:   PetscInt       i,j,Nvlocal;
573:   PetscReal      alpha,lambda;
574:   PetscScalar    *x,*f;
575:   VecScatter     scatter;
576:   Vec            localX = user->localX;
577: #if defined (UNUSED_VARIABLES)
578:   PetscScalar    ut,ub,ul,ur,u,*g,sc,uyy,uxx;
579:   PetscReal      hx,hy,hxdhy,hydhx;
580:   PetscReal      two = 2.0,one = 1.0;
581:   PetscInt       Nvglobal,Neglobal,row;
582:   PetscInt       *gloInd;

584:   Nvglobal = user->Nvglobal;
585:   Neglobal = user->Neglobal;
586:   gloInd   = user->gloInd;
587: #endif

589:   Nvlocal  = user->Nvlocal;
590:   lambda   = user->non_lin_param;
591:   alpha    = user->lin_param;
592:   scatter  = user->scatter;

594:   /* 
595:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
596:      described in the beginning of this code 
597:                                                                                    
598:      First scatter the distributed vector X into local vector localX (that includes
599:      values for ghost nodes. If we wish,we can put some other work between 
600:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
601:      computation.
602:  */
603:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
604:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);

606:   /*
607:      Get pointers to vector data
608:   */
609:   VecGetArray(localX,&x);
610:   VecGetArray(F,&f);

612:   /* 
613:     Now compute the f(x). As mentioned earlier, the computed Laplacian is just an 
614:     approximate one chosen for illustrative purpose only. Another point to notice 
615:     is that this is a local (completly parallel) calculation. In practical application 
616:     codes, function calculation time is a dominat portion of the overall execution time.
617:   */
618:   for (i=0; i < Nvlocal; i++) {
619:     f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i];
620:     for (j = 0; j < user->itot[i]; j++) {
621:       f[i] -= x[user->AdjM[i][j]];
622:     }
623:   }

625:   /*
626:      Restore vectors
627:   */
628:   VecRestoreArray(localX,&x);
629:   VecRestoreArray(F,&f);
630:   /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/

632:   return 0;
633: }

637: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
638: /*
639:    FormJacobian - Evaluates Jacobian matrix.

641:    Input Parameters:
642: .  snes - the SNES context
643: .  X - input vector
644: .  ptr - optional user-defined context, as set by SNESSetJacobian()

646:    Output Parameters:
647: .  A - Jacobian matrix
648: .  B - optionally different preconditioning matrix
649: .  flag - flag indicating matrix structure

651: */
652: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
653: {
654:   AppCtx *user = (AppCtx*)ptr;
655:   Mat     jac = *B;
656:   PetscInt     i,j,Nvlocal,col[50],ierr;
657:   PetscScalar  alpha,lambda,value[50];
658:   Vec          localX = user->localX;
659:   VecScatter scatter;
660:   PetscScalar  *x;
661: #if defined (UNUSED_VARIABLES)
662:   PetscScalar  two = 2.0,one = 1.0;
663:   PetscInt     row,Nvglobal,Neglobal;
664:   PetscInt     *gloInd;

666:   Nvglobal = user->Nvglobal;
667:   Neglobal = user->Neglobal;
668:   gloInd   = user->gloInd;
669: #endif
670: 
671:   /*printf("Entering into FormJacobian \n");*/
672:   Nvlocal  = user->Nvlocal;
673:   lambda   = user->non_lin_param;
674:   alpha    =  user->lin_param;
675:   scatter  = user->scatter;

677:   /* 
678:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
679:      described in the beginning of this code 
680:                                                                                    
681:      First scatter the distributed vector X into local vector localX (that includes
682:      values for ghost nodes. If we wish, we can put some other work between 
683:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
684:      computation.
685:   */
686:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
687:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
688: 
689:   /*
690:      Get pointer to vector data
691:   */
692:   VecGetArray(localX,&x);

694:   for (i=0; i < Nvlocal; i++) {
695:     col[0] = i;
696:     value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
697:     for (j = 0; j < user->itot[i]; j++) {
698:       col[j+1] = user->AdjM[i][j];
699:       value[j+1] = -1.0;
700:     }

702:   /* 
703:     Set the matrix values in the local ordering. Note that in order to use this
704:     feature we must call the routine MatSetLocalToGlobalMapping() after the 
705:     matrix has been created. 
706:   */
707:     MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);
708:   }

710:   /* 
711:      Assemble matrix, using the 2-step process:
712:        MatAssemblyBegin(), MatAssemblyEnd(). 
713:      Between these two calls, the pointer to vector data has been restored to
714:      demonstrate the use of overlapping communicationn with computation.
715:   */
716:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
717:   VecRestoreArray(localX,&x);
718:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

720:   /*
721:      Set flag to indicate that the Jacobian matrix retains an identical
722:      nonzero structure throughout all nonlinear iterations (although the
723:      values of the entries change). Thus, we can save some work in setting
724:      up the preconditioner (e.g., no need to redo symbolic factorization for
725:      ILU/ICC preconditioners).
726:       - If the nonzero structure of the matrix is different during
727:         successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
728:         must be used instead.  If you are unsure whether the matrix
729:         structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
730:       - Caution:  If you specify SAME_NONZERO_PATTERN, PETSc
731:         believes your assertion and does not check the structure
732:         of the matrix.  If you erroneously claim that the structure
733:         is the same when it actually is not, the new preconditioner
734:         will not function correctly.  Thus, use this optimization
735:         feature with caution!
736:   */
737:   *flag = SAME_NONZERO_PATTERN;

739:   /*
740:      Tell the matrix we will never add a new nonzero location to the
741:      matrix. If we do, it will generate an error.
742:   */
743:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
744:   /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
745:   return 0;
746: }
747: #else
748: #include <stdio.h>
749: int main(int argc,char **args)
750: {
751:   fprintf(stdout,"This example does not work for complex numbers.\n");
752:   return 0;
753: }
754: #endif