Actual source code: ex10.c


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

 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: #define MAX_ELEM      500  /* Maximum number of elements */
 44: #define MAX_VERT      100  /* Maximum number of vertices */
 45: #define MAX_VERT_ELEM   3  /* Vertices per element       */

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

 65: /*
 66:   User-defined routines
 67: */
 68: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
 69: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
 70: PetscErrorCode FormInitialGuess(AppCtx*,Vec);

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

106:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
107:      Initialize program
108:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

110:   PetscInitialize(&argc,&argv,"options.inf",help);
111:   MPI_Comm_rank(MPI_COMM_WORLD,&rank);
112:   MPI_Comm_size(MPI_COMM_WORLD,&size);

114:   /* The current input file options.inf is for 2 proc run only */

117:   /*
118:      Initialize problem parameters
119:   */
120:   user.Nvglobal = 16;      /*Global # of vertices  */
121:   user.Neglobal = 18;      /*Global # of elements  */

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

126:   user.non_lin_param = 0.06;

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

130:   user.lin_param = -1.0;

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

134:   user.Nvlocal = 0;
135:   user.Nelocal = 0;

137:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138:       Read the mesh and partitioning information
139:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

141:   /*
142:      Read the mesh and partitioning information from 'adj.in'.
143:      The file format is as follows.
144:      For each line the first entry is the processor rank where the
145:      current node belongs. The second entry is the number of
146:      neighbors of a node. The rest of the line is the adjacency
147:      list of a node. Currently this file is set up to work on two
148:      processors.

150:      This is not a very good example of reading input. In the future,
151:      we will put an example that shows the style that should be
152:      used in a real application, where partitioning will be done
153:      dynamically by calling partitioning routines (at present, we have
154:      a  ready interface to ParMeTiS).
155:    */
156:   fptr = fopen("adj.in","r");

159:   /*
160:      Each processor writes to the file output.<rank> where rank is the
161:      processor's rank.
162:   */
163:   sprintf(part_name,"output.%d",rank);
164:   fptr1 = fopen(part_name,"w");
166:   PetscMalloc1(user.Nvglobal,&user.gloInd);
167:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Rank is %d\n",rank);
168:   for (inode = 0; inode < user.Nvglobal; inode++) {
170:     sscanf(str,"%d",&dtmp);user.v2p[inode] = dtmp;
171:     if (user.v2p[inode] == rank) {
172:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %D belongs to processor %D\n",inode,user.v2p[inode]);

174:       user.gloInd[user.Nvlocal] = inode;
175:       sscanf(str,"%*d %d",&dtmp);
176:       nbrs = dtmp;
177:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"Number of neighbors for the vertex %D is %D\n",inode,nbrs);

179:       user.itot[user.Nvlocal] = nbrs;
180:       Nvneighborstotal       += nbrs;
181:       for (i = 0; i < user.itot[user.Nvlocal]; i++) {
182:         form[0]='\0';
183:         for (j=0; j < i+2; j++) {
184:           PetscStrlcat(form,"%*d ",sizeof(form));
185:         }
186:         PetscStrlcat(form,"%d",sizeof(form));

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

191:         PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[user.Nvlocal][i]);
192:       }
193:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
194:       user.Nvlocal++;
195:     }
196:   }
197:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"Total # of Local Vertices is %D \n",user.Nvlocal);

199:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200:      Create different orderings
201:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

203:   /*
204:     Create the local ordering list for vertices. First a list using the PETSc global
205:     ordering is created. Then we use the AO object to get the PETSc-to-application and
206:     application-to-PETSc mappings. Each vertex also gets a local index (stored in the
207:     locInd array).
208:   */
209:   MPI_Scan(&user.Nvlocal,&rstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
210:   rstart -= user.Nvlocal;
211:   PetscMalloc1(user.Nvlocal,&pordering);

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

215:   /*
216:     Create the AO object
217:   */
218:   AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
219:   PetscFree(pordering);

221:   /*
222:     Keep the global indices for later use
223:   */
224:   PetscMalloc1(user.Nvlocal,&user.locInd);
225:   PetscMalloc1(Nvneighborstotal,&tmp);

227:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228:     Demonstrate the use of AO functionality
229:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

235:     user.locInd[i] = user.gloInd[i];
236:   }
237:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
238:   jstart = 0;
239:   for (i=0; i < user.Nvlocal; i++) {
240:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.gloInd[i]);
241:     for (j=0; j < user.itot[i]; j++) {
242:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);

244:       tmp[j + jstart] = user.AdjM[i][j];
245:     }
246:     jstart += user.itot[i];
247:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
248:   }

250:   /*
251:     Now map the vlocal and neighbor lists to the PETSc ordering
252:   */
253:   AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);
254:   AOApplicationToPetsc(ao,Nvneighborstotal,tmp);
255:   AODestroy(&ao);

257:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"After AOApplicationToPetsc, local indices are : \n");
258:   for (i=0; i < user.Nvlocal; i++) {
259:     PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.locInd[i]);
260:   }
261:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");

263:   jstart = 0;
264:   for (i=0; i < user.Nvlocal; i++) {
265:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.locInd[i]);
266:     for (j=0; j < user.itot[i]; j++) {
267:       user.AdjM[i][j] = tmp[j+jstart];

269:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
270:     }
271:     jstart += user.itot[i];
272:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
273:   }

275:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
276:      Extract the ghost vertex information for each processor
277:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
278:   /*
279:    Next, we need to generate a list of vertices required for this processor
280:    and a local numbering scheme for all vertices required on this processor.
281:       vertices - integer array of all vertices needed on this processor in PETSc
282:                  global numbering; this list consists of first the "locally owned"
283:                  vertices followed by the ghost vertices.
284:       verticesmask - integer array that for each global vertex lists its local
285:                      vertex number (in vertices) + 1. If the global vertex is not
286:                      represented on this processor, then the corresponding
287:                      entry in verticesmask is zero

289:       Note: vertices and verticesmask are both Nvglobal in length; this may
290:     sound terribly non-scalable, but in fact is not so bad for a reasonable
291:     number of processors. Importantly, it allows us to use NO SEARCHING
292:     in setting up the data structures.
293:   */
294:   PetscMalloc1(user.Nvglobal,&vertices);
295:   PetscCalloc1(user.Nvglobal,&verticesmask);
296:   nvertices = 0;

298:   /*
299:     First load "owned vertices" into list
300:   */
301:   for (i=0; i < user.Nvlocal; i++) {
302:     vertices[nvertices++]        = user.locInd[i];
303:     verticesmask[user.locInd[i]] = nvertices;
304:   }

306:   /*
307:     Now load ghost vertices into list
308:   */
309:   for (i=0; i < user.Nvlocal; i++) {
310:     for (j=0; j < user.itot[i]; j++) {
311:       nb = user.AdjM[i][j];
312:       if (!verticesmask[nb]) {
313:         vertices[nvertices++] = nb;
314:         verticesmask[nb]      = nvertices;
315:       }
316:     }
317:   }

319:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
320:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n");
321:   for (i=0; i < nvertices; i++) {
322:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",vertices[i]);
323:   }
324:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");

326:   /*
327:      Map the vertices listed in the neighbors to the local numbering from
328:     the global ordering that they contained initially.
329:   */
330:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
331:   PetscFPrintf(PETSC_COMM_SELF,fptr1,"After mapping neighbors in the local contiguous ordering\n");
332:   for (i=0; i<user.Nvlocal; i++) {
333:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are :\n",i);
334:     for (j = 0; j < user.itot[i]; j++) {
335:       nb              = user.AdjM[i][j];
336:       user.AdjM[i][j] = verticesmask[nb] - 1;

338:       PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
339:     }
340:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
341:   }

343:   N = user.Nvglobal;

345:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
346:      Create vector and matrix data structures
347:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

349:   /*
350:     Create vector data structures
351:   */
352:   VecCreate(MPI_COMM_WORLD,&x);
353:   VecSetSizes(x,user.Nvlocal,N);
354:   VecSetFromOptions(x);
355:   VecDuplicate(x,&r);
356:   VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
357:   VecDuplicate(user.localX,&user.localF);

359:   /*
360:     Create the scatter between the global representation and the
361:     local representation
362:   */
363:   ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
364:   ISCreateBlock(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isglobal);
365:   VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);
366:   ISDestroy(&isglobal);
367:   ISDestroy(&islocal);

369:   /*
370:      Create matrix data structure; Just to keep the example simple, we have not done any
371:      preallocation of memory for the matrix. In real application code with big matrices,
372:      preallocation should always be done to expedite the matrix creation.
373:   */
374:   MatCreate(MPI_COMM_WORLD,&Jac);
375:   MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
376:   MatSetFromOptions(Jac);
377:   MatSetUp(Jac);

379:   /*
380:     The following routine allows us to set the matrix values in local ordering
381:   */
382:   ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isl2g);
383:   MatSetLocalToGlobalMapping(Jac,isl2g,isl2g);

385:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
386:      Create nonlinear solver context
387:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

389:   SNESCreate(MPI_COMM_WORLD,&snes);
390:   SNESSetType(snes,type);

392:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
393:     Set routines for function and Jacobian evaluation
394:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
395:   SNESSetFunction(snes,r,FormFunction,(void*)&user);

397:   PetscOptionsGetBool(NULL,NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,0);
398:   if (!fd_jacobian_coloring) {
399:     SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user);
400:   } else {  /* Use matfdcoloring */
401:     ISColoring   iscoloring;
402:     MatColoring  mc;

404:     /* Get the data structure of Jac */
405:     FormJacobian(snes,x,Jac,Jac,&user);
406:     /* Create coloring context */
407:     MatColoringCreate(Jac,&mc);
408:     MatColoringSetType(mc,MATCOLORINGSL);
409:     MatColoringSetFromOptions(mc);
410:     MatColoringApply(mc,&iscoloring);
411:     MatColoringDestroy(&mc);
412:     MatFDColoringCreate(Jac,iscoloring,&matfdcoloring);
413:     MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
414:     MatFDColoringSetFromOptions(matfdcoloring);
415:     MatFDColoringSetUp(Jac,iscoloring,matfdcoloring);
416:     /* MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD); */
417:     SNESSetJacobian(snes,Jac,Jac,SNESComputeJacobianDefaultColor,matfdcoloring);
418:     ISColoringDestroy(&iscoloring);
419:   }

421:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
422:      Customize nonlinear solver; set runtime options
423:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

425:   SNESSetFromOptions(snes);

427:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
428:      Evaluate initial guess; then solve nonlinear system
429:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

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

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

448:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
449:      Now solve the nonlinear system
450:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

452:   SNESSolve(snes,NULL,x);
453:   SNESGetIterationNumber(snes,&its);
454:   SNESGetNonlinearStepFailures(snes,&nfails);

456:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
457:     Print the output : solution vector and other information
458:      Each processor writes to the file output.<rank> where rank is the
459:      processor's rank.
460:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

462:   VecGetArray(x,&xx);
463:   for (inode = 0; inode < user.Nvlocal; inode++) {
464:     PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %D is %f \n",inode,xx[inode]);
465:   }
466:   VecRestoreArray(x,&xx);
467:   fclose(fptr1);
468:   PetscPrintf(MPI_COMM_WORLD,"number of SNES iterations = %D, ",its);
469:   PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %D\n",nfails);

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

495: /*
496:    FormInitialGuess - Forms initial approximation.

498:    Input Parameters:
499:    user - user-defined application context
500:    X - vector

502:    Output Parameter:
503:    X - vector
504:  */
505: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
506: {
507:   PetscInt    i,Nvlocal;
508:   PetscInt    *gloInd;
509:   PetscScalar *x;
510: #if defined(UNUSED_VARIABLES)
511:   PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc;
512:   PetscInt  Neglobal,Nvglobal,j,row;
513:   PetscReal alpha,lambda;

515:   Nvglobal = user->Nvglobal;
516:   Neglobal = user->Neglobal;
517:   lambda   = user->non_lin_param;
518:   alpha    = user->lin_param;
519: #endif

521:   Nvlocal = user->Nvlocal;
522:   gloInd  = user->gloInd;

524:   /*
525:      Get a pointer to vector data.
526:        - For default PETSc vectors, VecGetArray() returns a pointer to
527:          the data array.  Otherwise, the routine is implementation dependent.
528:        - You MUST call VecRestoreArray() when you no longer need access to
529:          the array.
530:   */
531:   VecGetArray(X,&x);

533:   /*
534:      Compute initial guess over the locally owned part of the grid
535:   */
536:   for (i=0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i];

538:   /*
539:      Restore vector
540:   */
541:   VecRestoreArray(X,&x);
542:   return 0;
543: }
544: /* --------------------  Evaluate Function F(x) --------------------- */
545: /*
546:    FormFunction - Evaluates nonlinear function, F(x).

548:    Input Parameters:
549: .  snes - the SNES context
550: .  X - input vector
551: .  ptr - optional user-defined context, as set by SNESSetFunction()

553:    Output Parameter:
554: .  F - function vector
555:  */
556: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
557: {
558:   AppCtx         *user = (AppCtx*)ptr;
559:   PetscInt       i,j,Nvlocal;
560:   PetscReal      alpha,lambda;
561:   PetscScalar    *x,*f;
562:   VecScatter     scatter;
563:   Vec            localX = user->localX;
564: #if defined(UNUSED_VARIABLES)
565:   PetscScalar ut,ub,ul,ur,u,*g,sc,uyy,uxx;
566:   PetscReal   hx,hy,hxdhy,hydhx;
567:   PetscReal   two = 2.0,one = 1.0;
568:   PetscInt    Nvglobal,Neglobal,row;
569:   PetscInt    *gloInd;

571:   Nvglobal = user->Nvglobal;
572:   Neglobal = user->Neglobal;
573:   gloInd   = user->gloInd;
574: #endif

576:   Nvlocal = user->Nvlocal;
577:   lambda  = user->non_lin_param;
578:   alpha   = user->lin_param;
579:   scatter = user->scatter;

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

585:      First scatter the distributed vector X into local vector localX (that includes
586:      values for ghost nodes. If we wish,we can put some other work between
587:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
588:      computation.
589:  */
590:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
591:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);

593:   /*
594:      Get pointers to vector data
595:   */
596:   VecGetArray(localX,&x);
597:   VecGetArray(F,&f);

599:   /*
600:     Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
601:     approximate one chosen for illustrative purpose only. Another point to notice
602:     is that this is a local (completly parallel) calculation. In practical application
603:     codes, function calculation time is a dominat portion of the overall execution time.
604:   */
605:   for (i=0; i < Nvlocal; i++) {
606:     f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i];
607:     for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]];
608:   }

610:   /*
611:      Restore vectors
612:   */
613:   VecRestoreArray(localX,&x);
614:   VecRestoreArray(F,&f);
615:   /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/

617:   return 0;
618: }

620: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
621: /*
622:    FormJacobian - Evaluates Jacobian matrix.

624:    Input Parameters:
625: .  snes - the SNES context
626: .  X - input vector
627: .  ptr - optional user-defined context, as set by SNESSetJacobian()

629:    Output Parameters:
630: .  A - Jacobian matrix
631: .  B - optionally different preconditioning matrix
632: .  flag - flag indicating matrix structure

634: */
635: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
636: {
637:   AppCtx      *user = (AppCtx*)ptr;
638:   PetscInt    i,j,Nvlocal,col[50];
639:   PetscScalar alpha,lambda,value[50];
640:   Vec         localX = user->localX;
641:   VecScatter  scatter;
642:   PetscScalar *x;
643: #if defined(UNUSED_VARIABLES)
644:   PetscScalar two = 2.0,one = 1.0;
645:   PetscInt    row,Nvglobal,Neglobal;
646:   PetscInt    *gloInd;

648:   Nvglobal = user->Nvglobal;
649:   Neglobal = user->Neglobal;
650:   gloInd   = user->gloInd;
651: #endif

653:   /*printf("Entering into FormJacobian \n");*/
654:   Nvlocal = user->Nvlocal;
655:   lambda  = user->non_lin_param;
656:   alpha   =  user->lin_param;
657:   scatter = user->scatter;

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

663:      First scatter the distributed vector X into local vector localX (that includes
664:      values for ghost nodes. If we wish, we can put some other work between
665:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
666:      computation.
667:   */
668:   VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
669:   VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);

671:   /*
672:      Get pointer to vector data
673:   */
674:   VecGetArray(localX,&x);

676:   for (i=0; i < Nvlocal; i++) {
677:     col[0]   = i;
678:     value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
679:     for (j = 0; j < user->itot[i]; j++) {
680:       col[j+1]   = user->AdjM[i][j];
681:       value[j+1] = -1.0;
682:     }

684:     /*
685:       Set the matrix values in the local ordering. Note that in order to use this
686:       feature we must call the routine MatSetLocalToGlobalMapping() after the
687:       matrix has been created.
688:     */
689:     MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);
690:   }

692:   /*
693:      Assemble matrix, using the 2-step process:
694:        MatAssemblyBegin(), MatAssemblyEnd().
695:      Between these two calls, the pointer to vector data has been restored to
696:      demonstrate the use of overlapping communicationn with computation.
697:   */
698:   MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
699:   VecRestoreArray(localX,&x);
700:   MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);

702:   /*
703:      Tell the matrix we will never add a new nonzero location to the
704:      matrix. If we do, it will generate an error.
705:   */
706:   MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
707:   /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
708:   return 0;
709: }

711: /*TEST

713:    build:
714:       requires: !complex

716:    test:
717:       nsize: 2
718:       args: -snes_monitor_short
719:       localrunfiles: options.inf adj.in

721:    test:
722:       suffix: 2
723:       nsize: 2
724:       args: -snes_monitor_short -fd_jacobian_coloring
725:       localrunfiles: options.inf adj.in

727: TEST*/