Actual source code: ex10.c
1: /*
2: Program usage: mpirun -np <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: The command line options include:\n\
19: -vert <Nv>, where Nv is the global number of nodes\n\
20: -elem <Ne>, where Ne is the global number of elements\n\
21: -nl_par <lambda>, where lambda is the multiplier for the non linear term (u*u) term\n\
22: -lin_par <alpha>, where alpha is the multiplier for the linear term (u) \n";
24: /*T
25: Concepts: SNES^unstructured grid
26: Concepts: AO^application to PETSc ordering or vice versa;
27: Concepts: VecScatter^using vector scatter operations;
28: Processors: n
29: T*/
31: /* ------------------------------------------------------------------------
33: PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian.
35: The Laplacian is approximated in the following way: each edge is given a weight
36: of one meaning that the diagonal term will have the weight equal to the degree
37: of a node. The off diagonal terms will get a weight of -1.
39: -----------------------------------------------------------------------*/
41: /*
42: Include petscao.h so that we can use AO (Application Ordering) object's services.
43: Include "petscsnes.h" so that we can use SNES solvers. Note that this
44: file automatically includes:
45: petsc.h - base PETSc routines petscvec.h - vectors
46: petscsys.h - system routines petscmat.h - matrices
47: petscis.h - index sets petscksp.h - Krylov subspace methods
48: petscviewer.h - viewers petscpc.h - preconditioners
49: petscksp.h - linear solvers
50: */
51: #include "petscao.h"
52: #include "petscsnes.h"
55: #define MAX_ELEM 500 /* Maximum number of elements */
56: #define MAX_VERT 100 /* Maximum number of vertices */
57: #define MAX_VERT_ELEM 3 /* Vertices per element */
59: /*
60: Application-defined context for problem specific data
61: */
62: typedef struct {
63: PetscInt Nvglobal,Nvlocal; /* global and local number of vertices */
64: PetscInt Neglobal,Nelocal; /* global and local number of vertices */
65: PetscInt AdjM[MAX_VERT][50]; /* adjacency list of a vertex */
66: PetscInt itot[MAX_VERT]; /* total number of neighbors for a vertex */
67: PetscInt icv[MAX_ELEM][MAX_VERT_ELEM]; /* vertices belonging to an element */
68: PetscInt v2p[MAX_VERT]; /* processor number for a vertex */
69: PetscInt *locInd,*gloInd; /* local and global orderings for a node */
70: Vec localX,localF; /* local solution (u) and f(u) vectors */
71: PetscReal non_lin_param; /* nonlinear parameter for the PDE */
72: PetscReal lin_param; /* linear parameter for the PDE */
73: VecScatter scatter; /* scatter context for the local and
74: distributed vectors */
75: } AppCtx;
77: /*
78: User-defined routines
79: */
80: PetscErrorCode FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*),
81: FormFunction(SNES,Vec,Vec,void*),
82: FormInitialGuess(AppCtx*,Vec);
86: int main(int argc,char **argv)
87: {
88: SNES snes; /* SNES context */
89: SNESType type = SNESLS; /* default nonlinear solution method */
90: Vec x,r; /* solution, residual vectors */
91: Mat Jac; /* Jacobian matrix */
92: AppCtx user; /* user-defined application context */
93: AO ao; /* Application Ordering object */
94: IS isglobal,islocal; /* global and local index sets */
95: PetscMPIInt rank,size; /* rank of a process, number of processors */
96: PetscInt rstart; /* starting index of PETSc ordering for a processor */
97: PetscInt nfails; /* number of unsuccessful Newton steps */
98: PetscInt bs = 1; /* block size for multicomponent systems */
99: PetscInt nvertices; /* number of local plus ghost nodes of a processor */
100: PetscInt *pordering; /* PETSc ordering */
101: PetscInt *vertices; /* list of all vertices (incl. ghost ones)
102: on a processor */
103: PetscInt *verticesmask,*svertices;
104: PetscInt *tmp;
105: PetscInt i,j,jstart,inode,nb,nbrs,Nvneighborstotal = 0;
106: PetscErrorCode ierr;
107: PetscInt its,N;
108: PetscScalar *xx;
109: char str[256],form[256],part_name[256];
110: FILE *fptr,*fptr1;
111: ISLocalToGlobalMapping isl2g;
112: int dtmp;
113: #if defined (UNUSED_VARIABLES)
114: PetscDraw draw; /* drawing context */
115: PetscScalar *ff,*gg;
116: PetscReal tiny = 1.0e-10,zero = 0.0,one = 1.0,big = 1.0e+10;
117: PetscInt *tmp1,*tmp2;
118: #endif
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(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 */
135: PetscOptionsGetInt(PETSC_NULL,"-vert",&user.Nvglobal,PETSC_NULL);
136: PetscOptionsGetInt(PETSC_NULL,"-elem",&user.Neglobal,PETSC_NULL);
137: user.non_lin_param = 0.06;
138: PetscOptionsGetReal(PETSC_NULL,"-nl_par",&user.non_lin_param,PETSC_NULL);
139: user.lin_param = -1.0;
140: PetscOptionsGetReal(PETSC_NULL,"-lin_par",&user.lin_param,PETSC_NULL);
141: user.Nvlocal = 0;
142: user.Nelocal = 0;
144: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
145: Read the mesh and partitioning information
146: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
147:
148: /*
149: Read the mesh and partitioning information from 'adj.in'.
150: The file format is as follows.
151: For each line the first entry is the processor rank where the
152: current node belongs. The second entry is the number of
153: neighbors of a node. The rest of the line is the adjacency
154: list of a node. Currently this file is set up to work on two
155: processors.
157: This is not a very good example of reading input. In the future,
158: we will put an example that shows the style that should be
159: used in a real application, where partitioning will be done
160: dynamically by calling partitioning routines (at present, we have
161: a ready interface to ParMeTiS).
162: */
163: fptr = fopen("adj.in","r");
164: if (!fptr) {
165: SETERRQ(0,"Could not open adj.in")
166: }
167:
168: /*
169: Each processor writes to the file output.<rank> where rank is the
170: processor's rank.
171: */
172: sprintf(part_name,"output.%d",rank);
173: fptr1 = fopen(part_name,"w");
174: if (!fptr1) {
175: SETERRQ(0,"Could no open output file");
176: }
177: PetscMalloc(user.Nvglobal*sizeof(PetscInt),&user.gloInd);
178: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Rank is %D\n",rank);
179: for (inode = 0; inode < user.Nvglobal; inode++) {
180: fgets(str,256,fptr);
181: sscanf(str,"%d",&dtmp);user.v2p[inode] = dtmp;
182: if (user.v2p[inode] == rank) {
183: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %D belongs to processor %D\n",inode,user.v2p[inode]);
184: user.gloInd[user.Nvlocal] = inode;
185: sscanf(str,"%*d %d",&dtmp);nbrs = dtmp;
186: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Number of neighbors for the vertex %D is %D\n",inode,nbrs);
187: user.itot[user.Nvlocal] = nbrs;
188: Nvneighborstotal += nbrs;
189: for (i = 0; i < user.itot[user.Nvlocal]; i++){
190: form[0]='\0';
191: for (j=0; j < i+2; j++){
192: PetscStrcat(form,"%*d ");
193: }
194: PetscStrcat(form,"%d");
195: sscanf(str,form,&dtmp); user.AdjM[user.Nvlocal][i] = dtmp;
196: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[user.Nvlocal][i]);
197: }
198: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
199: user.Nvlocal++;
200: }
201: }
202: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Total # of Local Vertices is %D \n",user.Nvlocal);
204: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
205: Create different orderings
206: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
208: /*
209: Create the local ordering list for vertices. First a list using the PETSc global
210: ordering is created. Then we use the AO object to get the PETSc-to-application and
211: application-to-PETSc mappings. Each vertex also gets a local index (stored in the
212: locInd array).
213: */
214: MPI_Scan(&user.Nvlocal,&rstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
215: rstart -= user.Nvlocal;
216: PetscMalloc(user.Nvlocal*sizeof(PetscInt),&pordering);
218: for (i=0; i < user.Nvlocal; i++) {
219: pordering[i] = rstart + i;
220: }
222: /*
223: Create the AO object
224: */
225: AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
226: PetscFree(pordering);
227:
228: /*
229: Keep the global indices for later use
230: */
231: PetscMalloc(user.Nvlocal*sizeof(PetscInt),&user.locInd);
232: PetscMalloc(Nvneighborstotal*sizeof(PetscInt),&tmp);
233:
234: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: Demonstrate the use of AO functionality
236: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
238: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Before AOApplicationToPetsc, local indices are : \n");
239: for (i=0; i < user.Nvlocal; i++) {
240: PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.gloInd[i]);
241: user.locInd[i] = user.gloInd[i];
242: }
243: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
244: jstart = 0;
245: for (i=0; i < user.Nvlocal; i++) {
246: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.gloInd[i]);
247: for (j=0; j < user.itot[i]; j++) {
248: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
249: tmp[j + jstart] = user.AdjM[i][j];
250: }
251: jstart += user.itot[i];
252: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
253: }
255: /*
256: Now map the vlocal and neighbor lists to the PETSc ordering
257: */
258: AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);
259: AOApplicationToPetsc(ao,Nvneighborstotal,tmp);
260: AODestroy(ao);
262: PetscFPrintf(PETSC_COMM_SELF,fptr1,"After AOApplicationToPetsc, local indices are : \n");
263: for (i=0; i < user.Nvlocal; i++) {
264: PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.locInd[i]);
265: }
266: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
268: jstart = 0;
269: for (i=0; i < user.Nvlocal; i++) {
270: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.locInd[i]);
271: for (j=0; j < user.itot[i]; j++) {
272: user.AdjM[i][j] = tmp[j+jstart];
273: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
274: }
275: jstart += user.itot[i];
276: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
277: }
279: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
280: Extract the ghost vertex information for each processor
281: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
282: /*
283: Next, we need to generate a list of vertices required for this processor
284: and a local numbering scheme for all vertices required on this processor.
285: vertices - integer array of all vertices needed on this processor in PETSc
286: global numbering; this list consists of first the "locally owned"
287: vertices followed by the ghost vertices.
288: verticesmask - integer array that for each global vertex lists its local
289: vertex number (in vertices) + 1. If the global vertex is not
290: represented on this processor, then the corresponding
291: entry in verticesmask is zero
292:
293: Note: vertices and verticesmask are both Nvglobal in length; this may
294: sound terribly non-scalable, but in fact is not so bad for a reasonable
295: number of processors. Importantly, it allows us to use NO SEARCHING
296: in setting up the data structures.
297: */
298: PetscMalloc(user.Nvglobal*sizeof(PetscInt),&vertices);
299: PetscMalloc(user.Nvglobal*sizeof(PetscInt),&verticesmask);
300: PetscMemzero(verticesmask,user.Nvglobal*sizeof(PetscInt));
301: nvertices = 0;
302:
303: /*
304: First load "owned vertices" into list
305: */
306: for (i=0; i < user.Nvlocal; i++) {
307: vertices[nvertices++] = user.locInd[i];
308: verticesmask[user.locInd[i]] = nvertices;
309: }
310:
311: /*
312: Now load ghost vertices into list
313: */
314: for (i=0; i < user.Nvlocal; i++) {
315: for (j=0; j < user.itot[i]; j++) {
316: nb = user.AdjM[i][j];
317: if (!verticesmask[nb]) {
318: vertices[nvertices++] = nb;
319: verticesmask[nb] = nvertices;
320: }
321: }
322: }
324: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
325: PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n");
326: for (i=0; i < nvertices; i++) {
327: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",vertices[i]);
328: }
329: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
330:
331: /*
332: Map the vertices listed in the neighbors to the local numbering from
333: the global ordering that they contained initially.
334: */
335: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
336: PetscFPrintf(PETSC_COMM_SELF,fptr1,"After mapping neighbors in the local contiguous ordering\n");
337: for (i=0; i<user.Nvlocal; i++) {
338: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are :\n",i);
339: for (j = 0; j < user.itot[i]; j++) {
340: nb = user.AdjM[i][j];
341: user.AdjM[i][j] = verticesmask[nb] - 1;
342: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
343: }
344: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
345: }
347: N = user.Nvglobal;
348:
349: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
350: Create vector and matrix data structures
351: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
353: /*
354: Create vector data structures
355: */
356: VecCreate(MPI_COMM_WORLD,&x);
357: VecSetSizes(x,user.Nvlocal,N);
358: VecSetFromOptions(x);
359: VecDuplicate(x,&r);
360: VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
361: VecDuplicate(user.localX,&user.localF);
363: /*
364: Create the scatter between the global representation and the
365: local representation
366: */
367: ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
368: PetscMalloc(nvertices*sizeof(PetscInt),&svertices);
369: for (i=0; i<nvertices; i++) svertices[i] = bs*vertices[i];
370: ISCreateBlock(MPI_COMM_SELF,bs,nvertices,svertices,&isglobal);
371: PetscFree(svertices);
372: VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);
374: /*
375: Create matrix data structure; Just to keep the example simple, we have not done any
376: preallocation of memory for the matrix. In real application code with big matrices,
377: preallocation should always be done to expedite the matrix creation.
378: */
379: MatCreate(MPI_COMM_WORLD,&Jac);
380: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
381: MatSetFromOptions(Jac);
383: /*
384: The following routine allows us to set the matrix values in local ordering
385: */
386: ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs*nvertices,vertices,&isl2g);
387: MatSetLocalToGlobalMapping(Jac,isl2g);
389: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
390: Create nonlinear solver context
391: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
393: SNESCreate(MPI_COMM_WORLD,&snes);
394: SNESSetType(snes,type);
396: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
397: Set routines for function and Jacobian evaluation
398: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
400: FormInitialGuess(&user,x);
401: SNESSetFunction(snes,r,FormFunction,(void*)&user);
402: SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user);
404: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
405: Customize nonlinear solver; set runtime options
406: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
408: SNESSetFromOptions(snes);
410: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
411: Evaluate initial guess; then solve nonlinear system
412: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
414: /*
415: Note: The user should initialize the vector, x, with the initial guess
416: for the nonlinear solver prior to calling SNESSolve(). In particular,
417: to employ an initial guess of zero, the user should explicitly set
418: this vector to zero by calling VecSet().
419: */
420: FormInitialGuess(&user,x);
422: /*
423: Print the initial guess
424: */
425: VecGetArray(x,&xx);
426: for (inode = 0; inode < user.Nvlocal; inode++)
427: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Initial Solution at node %D is %f \n",inode,xx[inode]);
428: VecRestoreArray(x,&xx);
430: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
431: Now solve the nonlinear system
432: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
434: SNESSolve(snes,PETSC_NULL,x);
435: SNESGetIterationNumber(snes,&its);
436: SNESGetNumberUnsuccessfulSteps(snes,&nfails);
437:
438: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
439: Print the output : solution vector and other information
440: Each processor writes to the file output.<rank> where rank is the
441: processor's rank.
442: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
444: VecGetArray(x,&xx);
445: for (inode = 0; inode < user.Nvlocal; inode++)
446: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %D is %f \n",inode,xx[inode]);
447: VecRestoreArray(x,&xx);
448: fclose(fptr1);
449: PetscPrintf(MPI_COMM_WORLD,"number of Newton iterations = %D, ",its);
450: PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %D\n",nfails);
452: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
453: Free work space. All PETSc objects should be destroyed when they
454: are no longer needed.
455: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
457: VecDestroy(x);
458: VecDestroy(r);
459: VecDestroy(user.localX);
460: VecDestroy(user.localF);
461: MatDestroy(Jac); SNESDestroy(snes);
462: /*PetscDrawDestroy(draw);*/
463: PetscFinalize();
465: return 0;
466: }
469: /* -------------------- Form initial approximation ----------------- */
471: /*
472: FormInitialGuess - Forms initial approximation.
474: Input Parameters:
475: user - user-defined application context
476: X - vector
478: Output Parameter:
479: X - vector
480: */
481: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
482: {
483: PetscInt i,Nvlocal,ierr;
484: PetscInt *gloInd;
485: PetscScalar *x;
486: #if defined (UNUSED_VARIABLES)
487: PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc;
488: PetscInt Neglobal,Nvglobal,j,row;
489: PetscReal alpha,lambda;
491: Nvglobal = user->Nvglobal;
492: Neglobal = user->Neglobal;
493: lambda = user->non_lin_param;
494: alpha = user->lin_param;
495: #endif
497: Nvlocal = user->Nvlocal;
498: gloInd = user->gloInd;
500: /*
501: Get a pointer to vector data.
502: - For default PETSc vectors, VecGetArray() returns a pointer to
503: the data array. Otherwise, the routine is implementation dependent.
504: - You MUST call VecRestoreArray() when you no longer need access to
505: the array.
506: */
507: VecGetArray(X,&x);
509: /*
510: Compute initial guess over the locally owned part of the grid
511: */
512: for (i=0; i < Nvlocal; i++) {
513: x[i] = (PetscReal)gloInd[i];
514: }
516: /*
517: Restore vector
518: */
519: VecRestoreArray(X,&x);
520: return 0;
521: }
524: /* -------------------- Evaluate Function F(x) --------------------- */
525: /*
526: FormFunction - Evaluates nonlinear function, F(x).
528: Input Parameters:
529: . snes - the SNES context
530: . X - input vector
531: . ptr - optional user-defined context, as set by SNESSetFunction()
533: Output Parameter:
534: . F - function vector
535: */
536: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
537: {
538: AppCtx *user = (AppCtx*)ptr;
540: PetscInt i,j,Nvlocal;
541: PetscReal alpha,lambda;
542: PetscScalar *x,*f;
543: VecScatter scatter;
544: Vec localX = user->localX;
545: #if defined (UNUSED_VARIABLES)
546: PetscScalar ut,ub,ul,ur,u,*g,sc,uyy,uxx;
547: PetscReal hx,hy,hxdhy,hydhx;
548: PetscReal two = 2.0,one = 1.0;
549: PetscInt Nvglobal,Neglobal,row;
550: PetscInt *gloInd;
552: Nvglobal = user->Nvglobal;
553: Neglobal = user->Neglobal;
554: gloInd = user->gloInd;
555: #endif
557: Nvlocal = user->Nvlocal;
558: lambda = user->non_lin_param;
559: alpha = user->lin_param;
560: scatter = user->scatter;
562: /*
563: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
564: described in the beginning of this code
565:
566: First scatter the distributed vector X into local vector localX (that includes
567: values for ghost nodes. If we wish,we can put some other work between
568: VecScatterBegin() and VecScatterEnd() to overlap the communication with
569: computation.
570: */
571: VecScatterBegin(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
572: VecScatterEnd(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
574: /*
575: Get pointers to vector data
576: */
577: VecGetArray(localX,&x);
578: VecGetArray(F,&f);
580: /*
581: Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
582: approximate one chosen for illustrative purpose only. Another point to notice
583: is that this is a local (completly parallel) calculation. In practical application
584: codes, function calculation time is a dominat portion of the overall execution time.
585: */
586: for (i=0; i < Nvlocal; i++) {
587: f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i];
588: for (j = 0; j < user->itot[i]; j++) {
589: f[i] -= x[user->AdjM[i][j]];
590: }
591: }
593: /*
594: Restore vectors
595: */
596: VecRestoreArray(localX,&x);
597: VecRestoreArray(F,&f);
598: /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/
600: return 0;
601: }
605: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
606: /*
607: FormJacobian - Evaluates Jacobian matrix.
609: Input Parameters:
610: . snes - the SNES context
611: . X - input vector
612: . ptr - optional user-defined context, as set by SNESSetJacobian()
614: Output Parameters:
615: . A - Jacobian matrix
616: . B - optionally different preconditioning matrix
617: . flag - flag indicating matrix structure
619: */
620: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat *J,Mat *B,MatStructure *flag,void *ptr)
621: {
622: AppCtx *user = (AppCtx*)ptr;
623: Mat jac = *B;
624: PetscInt i,j,Nvlocal,col[50],ierr;
625: PetscScalar alpha,lambda,value[50];
626: Vec localX = user->localX;
627: VecScatter scatter;
628: PetscScalar *x;
629: #if defined (UNUSED_VARIABLES)
630: PetscScalar two = 2.0,one = 1.0;
631: PetscInt row,Nvglobal,Neglobal;
632: PetscInt *gloInd;
634: Nvglobal = user->Nvglobal;
635: Neglobal = user->Neglobal;
636: gloInd = user->gloInd;
637: #endif
638:
639: /*printf("Entering into FormJacobian \n");*/
640: Nvlocal = user->Nvlocal;
641: lambda = user->non_lin_param;
642: alpha = user->lin_param;
643: scatter = user->scatter;
645: /*
646: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
647: described in the beginning of this code
648:
649: First scatter the distributed vector X into local vector localX (that includes
650: values for ghost nodes. If we wish, we can put some other work between
651: VecScatterBegin() and VecScatterEnd() to overlap the communication with
652: computation.
653: */
654: VecScatterBegin(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
655: VecScatterEnd(X,localX,INSERT_VALUES,SCATTER_FORWARD,scatter);
656:
657: /*
658: Get pointer to vector data
659: */
660: VecGetArray(localX,&x);
662: for (i=0; i < Nvlocal; i++) {
663: col[0] = i;
664: value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
665: for (j = 0; j < user->itot[i]; j++) {
666: col[j+1] = user->AdjM[i][j];
667: value[j+1] = -1.0;
668: }
670: /*
671: Set the matrix values in the local ordering. Note that in order to use this
672: feature we must call the routine MatSetLocalToGlobalMapping() after the
673: matrix has been created.
674: */
675: MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);
676: }
678: /*
679: Assemble matrix, using the 2-step process:
680: MatAssemblyBegin(), MatAssemblyEnd().
681: Between these two calls, the pointer to vector data has been restored to
682: demonstrate the use of overlapping communicationn with computation.
683: */
684: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
685: VecRestoreArray(localX,&x);
686: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
688: /*
689: Set flag to indicate that the Jacobian matrix retains an identical
690: nonzero structure throughout all nonlinear iterations (although the
691: values of the entries change). Thus, we can save some work in setting
692: up the preconditioner (e.g., no need to redo symbolic factorization for
693: ILU/ICC preconditioners).
694: - If the nonzero structure of the matrix is different during
695: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
696: must be used instead. If you are unsure whether the matrix
697: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
698: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
699: believes your assertion and does not check the structure
700: of the matrix. If you erroneously claim that the structure
701: is the same when it actually is not, the new preconditioner
702: will not function correctly. Thus, use this optimization
703: feature with caution!
704: */
705: *flag = SAME_NONZERO_PATTERN;
707: /*
708: Tell the matrix we will never add a new nonzero location to the
709: matrix. If we do, it will generate an error.
710: */
711: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR);
712: /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
713: return 0;
714: }
715: #else
716: #include <stdio.h>
717: int main(int argc,char **args)
718: {
719: fprintf(stdout,"This example does not work for complex numbers.\n");
720: return 0;
721: }
722: #endif