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*/