Actual source code: ex10.c

petsc-3.5.4 2015-05-23
Report Typos and Errors
  2: /*
  3:    Include "petscsnes.h" so that we can use SNES solvers.  Note that this
  4:    file automatically includes:
  5:      petscsys.h       - base PETSc routines   petscvec.h - vectors
  6:      petscmat.h - matrices
  7:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
  8:      petscviewer.h - viewers               petscpc.h  - preconditioners
  9:      petscksp.h   - linear solvers
 10: */
 11: #include <petscsnes.h>
 12: #include <petscao.h>

 14: #if !defined(PETSC_USE_COMPLEX)

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

 35: /*T
 36:    Concepts: SNES^unstructured grid
 37:    Concepts: AO^application to PETSc ordering or vice versa;
 38:    Concepts: VecScatter^using vector scatter operations;
 39:    Processors: n
 40: T*/

 42: /* ------------------------------------------------------------------------

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

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

 50:    -----------------------------------------------------------------------*/


 53: #define MAX_ELEM      500  /* Maximum number of elements */
 54: #define MAX_VERT      100  /* Maximum number of vertices */
 55: #define MAX_VERT_ELEM   3  /* Vertices per element       */

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

 75: /*
 76:   User-defined routines
 77: */
 78: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 79: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 80: PetscErrorCode FormInitialGuess(AppCtx*,Vec);

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

119:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
120:      Initialize program
121:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

123:   PetscInitialize(&argc,&argv,"options.inf",help);
124:   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
125:   MPI_Comm_size(MPI_COMM_WORLD,&size);

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

130:   /*
131:      Initialize problem parameters
132:   */
133:   user.Nvglobal = 16;      /*Global # of vertices  */
134:   user.Neglobal = 18;      /*Global # of elements  */

136:   PetscOptionsGetInt(NULL,"-vert",&user.Nvglobal,NULL);
137:   PetscOptionsGetInt(NULL,"-elem",&user.Neglobal,NULL);

139:   user.non_lin_param = 0.06;

141:   PetscOptionsGetReal(NULL,"-nl_par",&user.non_lin_param,NULL);

143:   user.lin_param = -1.0;

145:   PetscOptionsGetReal(NULL,"-lin_par",&user.lin_param,NULL);

147:   user.Nvlocal = 0;
148:   user.Nelocal = 0;

150:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
151:       Read the mesh and partitioning information
152:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

163:      This is not a very good example of reading input. In the future,
164:      we will put an example that shows the style that should be
165:      used in a real application, where partitioning will be done
166:      dynamically by calling partitioning routines (at present, we have
167:      a  ready interface to ParMeTiS).
168:    */
169:   fptr = fopen("adj.in","r");
170:   if (!fptr) SETERRQ(PETSC_COMM_SELF,0,"Could not open adj.in");

172:   /*
173:      Each processor writes to the file output.<rank> where rank is the
174:      processor's rank.
175:   */
176:   sprintf(part_name,"output.%d",rank);
177:   fptr1 = fopen(part_name,"w");
178:   if (!fptr1) SETERRQ(PETSC_COMM_SELF,0,"Could no open output file");
179:   PetscMalloc1(user.Nvglobal,&user.gloInd);
180:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Rank is %D\n",rank);
181:   for (inode = 0; inode < user.Nvglobal; inode++) {
182:     if (!fgets(str,256,fptr)) SETERRQ(PETSC_COMM_SELF,1,"fgets read failed");
183:     sscanf(str,"%d",&dtmp);user.v2p[inode] = dtmp;
184:     if (user.v2p[inode] == rank) {
185:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %D belongs to processor %D\n",inode,user.v2p[inode]);

187:       user.gloInd[user.Nvlocal] = inode;
188:       sscanf(str,"%*d %d",&dtmp);
189:       nbrs = dtmp;
190:       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");

201:         sscanf(str,form,&dtmp);
202:         user.AdjM[user.Nvlocal][i] = dtmp;

204:         PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[user.Nvlocal][i]);
205:       }
206:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
207:       user.Nvlocal++;
208:     }
209:   }
210:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Total # of Local Vertices is %D \n",user.Nvlocal);

212:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
213:      Create different orderings
214:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

226:   for (i=0; i < user.Nvlocal; i++) pordering[i] = rstart + i;

228:   /*
229:     Create the AO object
230:   */
231:   AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
232:   PetscFree(pordering);

234:   /*
235:     Keep the global indices for later use
236:   */
237:   PetscMalloc1(user.Nvlocal,&user.locInd);
238:   PetscMalloc1(Nvneighborstotal,&tmp);

240:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241:     Demonstrate the use of AO functionality
242:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

244:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Before AOApplicationToPetsc, local indices are : \n");
245:   for (i=0; i < user.Nvlocal; i++) {
246:     PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.gloInd[i]);

248:     user.locInd[i] = user.gloInd[i];
249:   }
250:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
251:   jstart = 0;
252:   for (i=0; i < user.Nvlocal; i++) {
253:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.gloInd[i]);
254:     for (j=0; j < user.itot[i]; j++) {
255:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);

257:       tmp[j + jstart] = user.AdjM[i][j];
258:     }
259:     jstart += user.itot[i];
260:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
261:   }

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

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

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

282:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
283:     }
284:     jstart += user.itot[i];
285:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
286:   }

288:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
289:      Extract the ghost vertex information for each processor
290:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
291:   /*
292:    Next, we need to generate a list of vertices required for this processor
293:    and a local numbering scheme for all vertices required on this processor.
294:       vertices - integer array of all vertices needed on this processor in PETSc
295:                  global numbering; this list consists of first the "locally owned"
296:                  vertices followed by the ghost vertices.
297:       verticesmask - integer array that for each global vertex lists its local
298:                      vertex number (in vertices) + 1. If the global vertex is not
299:                      represented on this processor, then the corresponding
300:                      entry in verticesmask is zero

302:       Note: vertices and verticesmask are both Nvglobal in length; this may
303:     sound terribly non-scalable, but in fact is not so bad for a reasonable
304:     number of processors. Importantly, it allows us to use NO SEARCHING
305:     in setting up the data structures.
306:   */
307:   PetscMalloc1(user.Nvglobal,&vertices);
308:   PetscMalloc1(user.Nvglobal,&verticesmask);
309:   PetscMemzero(verticesmask,user.Nvglobal*sizeof(PetscInt));
310:   nvertices = 0;

312:   /*
313:     First load "owned vertices" into list
314:   */
315:   for (i=0; i < user.Nvlocal; i++) {
316:     vertices[nvertices++]        = user.locInd[i];
317:     verticesmask[user.locInd[i]] = nvertices;
318:   }

320:   /*
321:     Now load ghost vertices into list
322:   */
323:   for (i=0; i < user.Nvlocal; i++) {
324:     for (j=0; j < user.itot[i]; j++) {
325:       nb = user.AdjM[i][j];
326:       if (!verticesmask[nb]) {
327:         vertices[nvertices++] = nb;
328:         verticesmask[nb]      = nvertices;
329:       }
330:     }
331:   }

333:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
334:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n");
335:   for (i=0; i < nvertices; i++) {
336:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",vertices[i]);
337:   }
338:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");

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

352:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
353:     }
354:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
355:   }

357:   N = user.Nvglobal;

359:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
360:      Create vector and matrix data structures
361:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

363:   /*
364:     Create vector data structures
365:   */
366:   VecCreate(MPI_COMM_WORLD,&x);
367:   VecSetSizes(x,user.Nvlocal,N);
368:   VecSetFromOptions(x);
369:   VecDuplicate(x,&r);
370:   VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
371:   VecDuplicate(user.localX,&user.localF);

373:   /*
374:     Create the scatter between the global representation and the
375:     local representation
376:   */
377:   ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
378:   ISCreateBlock(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isglobal);
379:   VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);
380:   ISDestroy(&isglobal);
381:   ISDestroy(&islocal);

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

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

399:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
400:      Create nonlinear solver context
401:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

403:   SNESCreate(MPI_COMM_WORLD,&snes);
404:   SNESSetType(snes,type);

406:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
407:     Set routines for function and Jacobian evaluation
408:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
409:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

411:   PetscOptionsGetBool(NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,0);
412:   if (!fd_jacobian_coloring) {
413:     SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user);
414:   } else {  /* Use matfdcoloring */
415:     ISColoring   iscoloring;
416:     MatColoring  mc;

418:     /* Get the data structure of Jac */
419:     FormJacobian(snes,x,Jac,Jac,&user);
420:     /* Create coloring context */
421:     MatColoringCreate(Jac,&mc);
422:     MatColoringSetType(mc,MATCOLORINGSL);
423:     MatColoringSetFromOptions(mc);
424:     MatColoringApply(mc,&iscoloring);
425:     MatColoringDestroy(&mc);
426:     MatFDColoringCreate(Jac,iscoloring,&matfdcoloring);
427:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
428:     MatFDColoringSetFromOptions(matfdcoloring);
429:     MatFDColoringSetUp(Jac,iscoloring,matfdcoloring);
430:     /* MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD); */
431:     SNESSetJacobian(snes,Jac,Jac,SNESComputeJacobianDefaultColor,matfdcoloring);
432:     ISColoringDestroy(&iscoloring);
433:   }

435:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
436:      Customize nonlinear solver; set runtime options
437:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

439:   SNESSetFromOptions(snes);

441:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
442:      Evaluate initial guess; then solve nonlinear system
443:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

445:   /*
446:      Note: The user should initialize the vector, x, with the initial guess
447:      for the nonlinear solver prior to calling SNESSolve().  In particular,
448:      to employ an initial guess of zero, the user should explicitly set
449:      this vector to zero by calling VecSet().
450:   */
451:   FormInitialGuess(&user,x);

453:   /*
454:     Print the initial guess
455:   */
456:   VecGetArray(x,&xx);
457:   for (inode = 0; inode < user.Nvlocal; inode++) {
458:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Initial Solution at node %D is %f \n",inode,xx[inode]);
459:   }
460:   VecRestoreArray(x,&xx);

462:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
463:      Now solve the nonlinear system
464:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

466:   SNESSolve(snes,NULL,x);
467:   SNESGetIterationNumber(snes,&its);
468:   SNESGetNonlinearStepFailures(snes,&nfails);

470:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
471:     Print the output : solution vector and other information
472:      Each processor writes to the file output.<rank> where rank is the
473:      processor's rank.
474:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

476:   VecGetArray(x,&xx);
477:   for (inode = 0; inode < user.Nvlocal; inode++) {
478:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %D is %f \n",inode,xx[inode]);
479:   }
480:   VecRestoreArray(x,&xx);
481:   fclose(fptr1);
482:   PetscPrintf(MPI_COMM_WORLD,"number of SNES iterations = %D, ",its);
483:   PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %D\n",nfails);

485:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
486:      Free work space.  All PETSc objects should be destroyed when they
487:      are no longer needed.
488:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
489:   PetscFree(user.gloInd);
490:   PetscFree(user.locInd);
491:   PetscFree(vertices);
492:   PetscFree(verticesmask);
493:   PetscFree(tmp);
494:   VecScatterDestroy(&user.scatter);
495:   ISLocalToGlobalMappingDestroy(&isl2g);
496:   VecDestroy(&x);
497:   VecDestroy(&r);
498:   VecDestroy(&user.localX);
499:   VecDestroy(&user.localF);
500:   MatDestroy(&Jac);  SNESDestroy(&snes);
501:   /*PetscDrawDestroy(draw);*/
502:   if (fd_jacobian_coloring) {
503:     MatFDColoringDestroy(&matfdcoloring);
504:   }
505:   PetscFinalize();

507:   return 0;
508: }
511: /* --------------------  Form initial approximation ----------------- */

513: /*
514:    FormInitialGuess - Forms initial approximation.

516:    Input Parameters:
517:    user - user-defined application context
518:    X - vector

520:    Output Parameter:
521:    X - vector
522:  */
523: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
524: {
525:   PetscInt    i,Nvlocal,ierr;
526:   PetscInt    *gloInd;
527:   PetscScalar *x;
528: #if defined(UNUSED_VARIABLES)
529:   PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc;
530:   PetscInt  Neglobal,Nvglobal,j,row;
531:   PetscReal alpha,lambda;

533:   Nvglobal = user->Nvglobal;
534:   Neglobal = user->Neglobal;
535:   lambda   = user->non_lin_param;
536:   alpha    = user->lin_param;
537: #endif

539:   Nvlocal = user->Nvlocal;
540:   gloInd  = user->gloInd;

542:   /*
543:      Get a pointer to vector data.
544:        - For default PETSc vectors, VecGetArray() returns a pointer to
545:          the data array.  Otherwise, the routine is implementation dependent.
546:        - You MUST call VecRestoreArray() when you no longer need access to
547:          the array.
548:   */
549:   VecGetArray(X,&x);

551:   /*
552:      Compute initial guess over the locally owned part of the grid
553:   */
554:   for (i=0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i];

556:   /*
557:      Restore vector
558:   */
559:   VecRestoreArray(X,&x);
560:   return 0;
561: }
564: /* --------------------  Evaluate Function F(x) --------------------- */
565: /*
566:    FormFunction - Evaluates nonlinear function, F(x).

568:    Input Parameters:
569: .  snes - the SNES context
570: .  X - input vector
571: .  ptr - optional user-defined context, as set by SNESSetFunction()

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

592:   Nvglobal = user->Nvglobal;
593:   Neglobal = user->Neglobal;
594:   gloInd   = user->gloInd;
595: #endif

597:   Nvlocal = user->Nvlocal;
598:   lambda  = user->non_lin_param;
599:   alpha   = user->lin_param;
600:   scatter = user->scatter;

602:   /*
603:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
604:      described in the beginning of this code

606:      First scatter the distributed vector X into local vector localX (that includes
607:      values for ghost nodes. If we wish,we can put some other work between
608:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
609:      computation.
610:  */
611:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
612:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);

614:   /*
615:      Get pointers to vector data
616:   */
617:   VecGetArray(localX,&x);
618:   VecGetArray(F,&f);

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

631:   /*
632:      Restore vectors
633:   */
634:   VecRestoreArray(localX,&x);
635:   VecRestoreArray(F,&f);
636:   /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/

638:   return 0;
639: }

643: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
644: /*
645:    FormJacobian - Evaluates Jacobian matrix.

647:    Input Parameters:
648: .  snes - the SNES context
649: .  X - input vector
650: .  ptr - optional user-defined context, as set by SNESSetJacobian()

652:    Output Parameters:
653: .  A - Jacobian matrix
654: .  B - optionally different preconditioning matrix
655: .  flag - flag indicating matrix structure

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

671:   Nvglobal = user->Nvglobal;
672:   Neglobal = user->Neglobal;
673:   gloInd   = user->gloInd;
674: #endif

676:   /*printf("Entering into FormJacobian \n");*/
677:   Nvlocal = user->Nvlocal;
678:   lambda  = user->non_lin_param;
679:   alpha   =  user->lin_param;
680:   scatter = user->scatter;

682:   /*
683:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
684:      described in the beginning of this code

686:      First scatter the distributed vector X into local vector localX (that includes
687:      values for ghost nodes. If we wish, we can put some other work between
688:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
689:      computation.
690:   */
691:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
692:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);

694:   /*
695:      Get pointer to vector data
696:   */
697:   VecGetArray(localX,&x);

699:   for (i=0; i < Nvlocal; i++) {
700:     col[0]   = i;
701:     value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
702:     for (j = 0; j < user->itot[i]; j++) {
703:       col[j+1]   = user->AdjM[i][j];
704:       value[j+1] = -1.0;
705:     }

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

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

725:   /*
726:      Tell the matrix we will never add a new nonzero location to the
727:      matrix. If we do, it will generate an error.
728:   */
729:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
730:   /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
731:   return 0;
732: }
733: #else

735: int main(int argc,char **args)
736: {
737:   fprintf(stdout,"This example does not work for complex numbers.\n");
738:   return 0;
739: }
740: #endif