Actual source code: ex10.c
petsc-3.8.4 2018-03-24
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/examples/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: /*T
34: Concepts: SNES^unstructured grid
35: Concepts: AO^application to PETSc ordering or vice versa;
36: Concepts: VecScatter^using vector scatter operations;
37: Processors: n
38: T*/
40: /* ------------------------------------------------------------------------
42: PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian.
44: The Laplacian is approximated in the following way: each edge is given a weight
45: of one meaning that the diagonal term will have the weight equal to the degree
46: of a node. The off diagonal terms will get a weight of -1.
48: -----------------------------------------------------------------------*/
51: #define MAX_ELEM 500 /* Maximum number of elements */
52: #define MAX_VERT 100 /* Maximum number of vertices */
53: #define MAX_VERT_ELEM 3 /* Vertices per element */
55: /*
56: Application-defined context for problem specific data
57: */
58: typedef struct {
59: PetscInt Nvglobal,Nvlocal; /* global and local number of vertices */
60: PetscInt Neglobal,Nelocal; /* global and local number of vertices */
61: PetscInt AdjM[MAX_VERT][50]; /* adjacency list of a vertex */
62: PetscInt itot[MAX_VERT]; /* total number of neighbors for a vertex */
63: PetscInt icv[MAX_ELEM][MAX_VERT_ELEM]; /* vertices belonging to an element */
64: PetscInt v2p[MAX_VERT]; /* processor number for a vertex */
65: PetscInt *locInd,*gloInd; /* local and global orderings for a node */
66: Vec localX,localF; /* local solution (u) and f(u) vectors */
67: PetscReal non_lin_param; /* nonlinear parameter for the PDE */
68: PetscReal lin_param; /* linear parameter for the PDE */
69: VecScatter scatter; /* scatter context for the local and
70: distributed vectors */
71: } AppCtx;
73: /*
74: User-defined routines
75: */
76: PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*);
77: PetscErrorCode FormFunction(SNES,Vec,Vec,void*);
78: PetscErrorCode FormInitialGuess(AppCtx*,Vec);
80: int main(int argc,char **argv)
81: {
82: SNES snes; /* SNES context */
83: SNESType type = SNESNEWTONLS; /* default nonlinear solution method */
84: Vec x,r; /* solution, residual vectors */
85: Mat Jac; /* Jacobian matrix */
86: AppCtx user; /* user-defined application context */
87: AO ao; /* Application Ordering object */
88: IS isglobal,islocal; /* global and local index sets */
89: PetscMPIInt rank,size; /* rank of a process, number of processors */
90: PetscInt rstart; /* starting index of PETSc ordering for a processor */
91: PetscInt nfails; /* number of unsuccessful Newton steps */
92: PetscInt bs = 1; /* block size for multicomponent systems */
93: PetscInt nvertices; /* number of local plus ghost nodes of a processor */
94: PetscInt *pordering; /* PETSc ordering */
95: PetscInt *vertices; /* list of all vertices (incl. ghost ones) on a processor */
96: PetscInt *verticesmask;
97: PetscInt *tmp;
98: PetscInt i,j,jstart,inode,nb,nbrs,Nvneighborstotal = 0;
99: PetscErrorCode ierr;
100: PetscInt its,N;
101: PetscScalar *xx;
102: char str[256],form[256],part_name[256];
103: FILE *fptr,*fptr1;
104: ISLocalToGlobalMapping isl2g;
105: int dtmp;
106: #if defined(UNUSED_VARIABLES)
107: PetscDraw draw; /* drawing context */
108: PetscScalar *ff,*gg;
109: PetscReal tiny = 1.0e-10,zero = 0.0,one = 1.0,big = 1.0e+10;
110: PetscInt *tmp1,*tmp2;
111: #endif
112: MatFDColoring matfdcoloring = 0;
113: PetscBool fd_jacobian_coloring = PETSC_FALSE;
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Initialize program
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: PetscInitialize(&argc,&argv,"options.inf",help);if (ierr) return ierr;
120: MPI_Comm_rank(MPI_COMM_WORLD,&rank);
121: MPI_Comm_size(MPI_COMM_WORLD,&size);
123: /* The current input file options.inf is for 2 proc run only */
124: if (size != 2) SETERRQ(PETSC_COMM_SELF,1,"This Example currently runs on 2 procs only.");
126: /*
127: Initialize problem parameters
128: */
129: user.Nvglobal = 16; /*Global # of vertices */
130: user.Neglobal = 18; /*Global # of elements */
132: PetscOptionsGetInt(NULL,NULL,"-vert",&user.Nvglobal,NULL);
133: PetscOptionsGetInt(NULL,NULL,"-elem",&user.Neglobal,NULL);
135: user.non_lin_param = 0.06;
137: PetscOptionsGetReal(NULL,NULL,"-nl_par",&user.non_lin_param,NULL);
139: user.lin_param = -1.0;
141: PetscOptionsGetReal(NULL,NULL,"-lin_par",&user.lin_param,NULL);
143: user.Nvlocal = 0;
144: user.Nelocal = 0;
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Read the mesh and partitioning information
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
150: /*
151: Read the mesh and partitioning information from 'adj.in'.
152: The file format is as follows.
153: For each line the first entry is the processor rank where the
154: current node belongs. The second entry is the number of
155: neighbors of a node. The rest of the line is the adjacency
156: list of a node. Currently this file is set up to work on two
157: processors.
159: This is not a very good example of reading input. In the future,
160: we will put an example that shows the style that should be
161: used in a real application, where partitioning will be done
162: dynamically by calling partitioning routines (at present, we have
163: a ready interface to ParMeTiS).
164: */
165: fptr = fopen("adj.in","r");
166: if (!fptr) SETERRQ(PETSC_COMM_SELF,0,"Could not open adj.in");
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) SETERRQ(PETSC_COMM_SELF,0,"Could no open output file");
175: PetscMalloc1(user.Nvglobal,&user.gloInd);
176: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Rank is %D\n",rank);
177: for (inode = 0; inode < user.Nvglobal; inode++) {
178: if (!fgets(str,256,fptr)) SETERRQ(PETSC_COMM_SELF,1,"fgets read failed");
179: sscanf(str,"%d",&dtmp);user.v2p[inode] = dtmp;
180: if (user.v2p[inode] == rank) {
181: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %D belongs to processor %D\n",inode,user.v2p[inode]);
183: user.gloInd[user.Nvlocal] = inode;
184: sscanf(str,"%*d %d",&dtmp);
185: nbrs = dtmp;
186: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Number of neighbors for the vertex %D is %D\n",inode,nbrs);
188: user.itot[user.Nvlocal] = nbrs;
189: Nvneighborstotal += nbrs;
190: for (i = 0; i < user.itot[user.Nvlocal]; i++) {
191: form[0]='\0';
192: for (j=0; j < i+2; j++) {
193: PetscStrcat(form,"%*d ");
194: }
195: PetscStrcat(form,"%d");
197: sscanf(str,form,&dtmp);
198: user.AdjM[user.Nvlocal][i] = dtmp;
200: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[user.Nvlocal][i]);
201: }
202: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
203: user.Nvlocal++;
204: }
205: }
206: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Total # of Local Vertices is %D \n",user.Nvlocal);
208: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
209: Create different orderings
210: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
212: /*
213: Create the local ordering list for vertices. First a list using the PETSc global
214: ordering is created. Then we use the AO object to get the PETSc-to-application and
215: application-to-PETSc mappings. Each vertex also gets a local index (stored in the
216: locInd array).
217: */
218: MPI_Scan(&user.Nvlocal,&rstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);
219: rstart -= user.Nvlocal;
220: PetscMalloc1(user.Nvlocal,&pordering);
222: for (i=0; i < user.Nvlocal; i++) pordering[i] = rstart + i;
224: /*
225: Create the AO object
226: */
227: AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);
228: PetscFree(pordering);
230: /*
231: Keep the global indices for later use
232: */
233: PetscMalloc1(user.Nvlocal,&user.locInd);
234: PetscMalloc1(Nvneighborstotal,&tmp);
236: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: Demonstrate the use of AO functionality
238: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
240: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Before AOApplicationToPetsc, local indices are : \n");
241: for (i=0; i < user.Nvlocal; i++) {
242: PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.gloInd[i]);
244: user.locInd[i] = user.gloInd[i];
245: }
246: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
247: jstart = 0;
248: for (i=0; i < user.Nvlocal; i++) {
249: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.gloInd[i]);
250: for (j=0; j < user.itot[i]; j++) {
251: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
253: tmp[j + jstart] = user.AdjM[i][j];
254: }
255: jstart += user.itot[i];
256: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
257: }
259: /*
260: Now map the vlocal and neighbor lists to the PETSc ordering
261: */
262: AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);
263: AOApplicationToPetsc(ao,Nvneighborstotal,tmp);
264: AODestroy(&ao);
266: PetscFPrintf(PETSC_COMM_SELF,fptr1,"After AOApplicationToPetsc, local indices are : \n");
267: for (i=0; i < user.Nvlocal; i++) {
268: PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.locInd[i]);
269: }
270: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
272: jstart = 0;
273: for (i=0; i < user.Nvlocal; i++) {
274: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.locInd[i]);
275: for (j=0; j < user.itot[i]; j++) {
276: user.AdjM[i][j] = tmp[j+jstart];
278: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
279: }
280: jstart += user.itot[i];
281: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
282: }
284: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
285: Extract the ghost vertex information for each processor
286: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
287: /*
288: Next, we need to generate a list of vertices required for this processor
289: and a local numbering scheme for all vertices required on this processor.
290: vertices - integer array of all vertices needed on this processor in PETSc
291: global numbering; this list consists of first the "locally owned"
292: vertices followed by the ghost vertices.
293: verticesmask - integer array that for each global vertex lists its local
294: vertex number (in vertices) + 1. If the global vertex is not
295: represented on this processor, then the corresponding
296: entry in verticesmask is zero
298: Note: vertices and verticesmask are both Nvglobal in length; this may
299: sound terribly non-scalable, but in fact is not so bad for a reasonable
300: number of processors. Importantly, it allows us to use NO SEARCHING
301: in setting up the data structures.
302: */
303: PetscMalloc1(user.Nvglobal,&vertices);
304: PetscMalloc1(user.Nvglobal,&verticesmask);
305: PetscMemzero(verticesmask,user.Nvglobal*sizeof(PetscInt));
306: nvertices = 0;
308: /*
309: First load "owned vertices" into list
310: */
311: for (i=0; i < user.Nvlocal; i++) {
312: vertices[nvertices++] = user.locInd[i];
313: verticesmask[user.locInd[i]] = nvertices;
314: }
316: /*
317: Now load ghost vertices into list
318: */
319: for (i=0; i < user.Nvlocal; i++) {
320: for (j=0; j < user.itot[i]; j++) {
321: nb = user.AdjM[i][j];
322: if (!verticesmask[nb]) {
323: vertices[nvertices++] = nb;
324: verticesmask[nb] = nvertices;
325: }
326: }
327: }
329: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
330: PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n");
331: for (i=0; i < nvertices; i++) {
332: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",vertices[i]);
333: }
334: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
336: /*
337: Map the vertices listed in the neighbors to the local numbering from
338: the global ordering that they contained initially.
339: */
340: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
341: PetscFPrintf(PETSC_COMM_SELF,fptr1,"After mapping neighbors in the local contiguous ordering\n");
342: for (i=0; i<user.Nvlocal; i++) {
343: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are :\n",i);
344: for (j = 0; j < user.itot[i]; j++) {
345: nb = user.AdjM[i][j];
346: user.AdjM[i][j] = verticesmask[nb] - 1;
348: PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);
349: }
350: PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");
351: }
353: N = user.Nvglobal;
355: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
356: Create vector and matrix data structures
357: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
359: /*
360: Create vector data structures
361: */
362: VecCreate(MPI_COMM_WORLD,&x);
363: VecSetSizes(x,user.Nvlocal,N);
364: VecSetFromOptions(x);
365: VecDuplicate(x,&r);
366: VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);
367: VecDuplicate(user.localX,&user.localF);
369: /*
370: Create the scatter between the global representation and the
371: local representation
372: */
373: ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);
374: ISCreateBlock(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isglobal);
375: VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);
376: ISDestroy(&isglobal);
377: ISDestroy(&islocal);
379: /*
380: Create matrix data structure; Just to keep the example simple, we have not done any
381: preallocation of memory for the matrix. In real application code with big matrices,
382: preallocation should always be done to expedite the matrix creation.
383: */
384: MatCreate(MPI_COMM_WORLD,&Jac);
385: MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);
386: MatSetFromOptions(Jac);
387: MatSetUp(Jac);
389: /*
390: The following routine allows us to set the matrix values in local ordering
391: */
392: ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isl2g);
393: MatSetLocalToGlobalMapping(Jac,isl2g,isl2g);
395: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
396: Create nonlinear solver context
397: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
399: SNESCreate(MPI_COMM_WORLD,&snes);
400: SNESSetType(snes,type);
402: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
403: Set routines for function and Jacobian evaluation
404: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
405: SNESSetFunction(snes,r,FormFunction,(void*)&user);
407: PetscOptionsGetBool(NULL,NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,0);
408: if (!fd_jacobian_coloring) {
409: SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user);
410: } else { /* Use matfdcoloring */
411: ISColoring iscoloring;
412: MatColoring mc;
414: /* Get the data structure of Jac */
415: FormJacobian(snes,x,Jac,Jac,&user);
416: /* Create coloring context */
417: MatColoringCreate(Jac,&mc);
418: MatColoringSetType(mc,MATCOLORINGSL);
419: MatColoringSetFromOptions(mc);
420: MatColoringApply(mc,&iscoloring);
421: MatColoringDestroy(&mc);
422: MatFDColoringCreate(Jac,iscoloring,&matfdcoloring);
423: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);
424: MatFDColoringSetFromOptions(matfdcoloring);
425: MatFDColoringSetUp(Jac,iscoloring,matfdcoloring);
426: /* MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD); */
427: SNESSetJacobian(snes,Jac,Jac,SNESComputeJacobianDefaultColor,matfdcoloring);
428: ISColoringDestroy(&iscoloring);
429: }
431: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
432: Customize nonlinear solver; set runtime options
433: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
435: SNESSetFromOptions(snes);
437: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
438: Evaluate initial guess; then solve nonlinear system
439: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
441: /*
442: Note: The user should initialize the vector, x, with the initial guess
443: for the nonlinear solver prior to calling SNESSolve(). In particular,
444: to employ an initial guess of zero, the user should explicitly set
445: this vector to zero by calling VecSet().
446: */
447: FormInitialGuess(&user,x);
449: /*
450: Print the initial guess
451: */
452: VecGetArray(x,&xx);
453: for (inode = 0; inode < user.Nvlocal; inode++) {
454: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Initial Solution at node %D is %f \n",inode,xx[inode]);
455: }
456: VecRestoreArray(x,&xx);
458: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
459: Now solve the nonlinear system
460: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
462: SNESSolve(snes,NULL,x);
463: SNESGetIterationNumber(snes,&its);
464: SNESGetNonlinearStepFailures(snes,&nfails);
466: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
467: Print the output : solution vector and other information
468: Each processor writes to the file output.<rank> where rank is the
469: processor's rank.
470: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
472: VecGetArray(x,&xx);
473: for (inode = 0; inode < user.Nvlocal; inode++) {
474: PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %D is %f \n",inode,xx[inode]);
475: }
476: VecRestoreArray(x,&xx);
477: fclose(fptr1);
478: PetscPrintf(MPI_COMM_WORLD,"number of SNES iterations = %D, ",its);
479: PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %D\n",nfails);
481: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
482: Free work space. All PETSc objects should be destroyed when they
483: are no longer needed.
484: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
485: PetscFree(user.gloInd);
486: PetscFree(user.locInd);
487: PetscFree(vertices);
488: PetscFree(verticesmask);
489: PetscFree(tmp);
490: VecScatterDestroy(&user.scatter);
491: ISLocalToGlobalMappingDestroy(&isl2g);
492: VecDestroy(&x);
493: VecDestroy(&r);
494: VecDestroy(&user.localX);
495: VecDestroy(&user.localF);
496: MatDestroy(&Jac); SNESDestroy(&snes);
497: /*PetscDrawDestroy(draw);*/
498: if (fd_jacobian_coloring) {
499: MatFDColoringDestroy(&matfdcoloring);
500: }
501: PetscFinalize();
502: return ierr;
503: }
504: /* -------------------- Form initial approximation ----------------- */
506: /*
507: FormInitialGuess - Forms initial approximation.
509: Input Parameters:
510: user - user-defined application context
511: X - vector
513: Output Parameter:
514: X - vector
515: */
516: PetscErrorCode FormInitialGuess(AppCtx *user,Vec X)
517: {
518: PetscInt i,Nvlocal,ierr;
519: PetscInt *gloInd;
520: PetscScalar *x;
521: #if defined(UNUSED_VARIABLES)
522: PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc;
523: PetscInt Neglobal,Nvglobal,j,row;
524: PetscReal alpha,lambda;
526: Nvglobal = user->Nvglobal;
527: Neglobal = user->Neglobal;
528: lambda = user->non_lin_param;
529: alpha = user->lin_param;
530: #endif
532: Nvlocal = user->Nvlocal;
533: gloInd = user->gloInd;
535: /*
536: Get a pointer to vector data.
537: - For default PETSc vectors, VecGetArray() returns a pointer to
538: the data array. Otherwise, the routine is implementation dependent.
539: - You MUST call VecRestoreArray() when you no longer need access to
540: the array.
541: */
542: VecGetArray(X,&x);
544: /*
545: Compute initial guess over the locally owned part of the grid
546: */
547: for (i=0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i];
549: /*
550: Restore vector
551: */
552: VecRestoreArray(X,&x);
553: return 0;
554: }
555: /* -------------------- Evaluate Function F(x) --------------------- */
556: /*
557: FormFunction - Evaluates nonlinear function, F(x).
559: Input Parameters:
560: . snes - the SNES context
561: . X - input vector
562: . ptr - optional user-defined context, as set by SNESSetFunction()
564: Output Parameter:
565: . F - function vector
566: */
567: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr)
568: {
569: AppCtx *user = (AppCtx*)ptr;
571: PetscInt i,j,Nvlocal;
572: PetscReal alpha,lambda;
573: PetscScalar *x,*f;
574: VecScatter scatter;
575: Vec localX = user->localX;
576: #if defined(UNUSED_VARIABLES)
577: PetscScalar ut,ub,ul,ur,u,*g,sc,uyy,uxx;
578: PetscReal hx,hy,hxdhy,hydhx;
579: PetscReal two = 2.0,one = 1.0;
580: PetscInt Nvglobal,Neglobal,row;
581: PetscInt *gloInd;
583: Nvglobal = user->Nvglobal;
584: Neglobal = user->Neglobal;
585: gloInd = user->gloInd;
586: #endif
588: Nvlocal = user->Nvlocal;
589: lambda = user->non_lin_param;
590: alpha = user->lin_param;
591: scatter = user->scatter;
593: /*
594: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
595: described in the beginning of this code
597: First scatter the distributed vector X into local vector localX (that includes
598: values for ghost nodes. If we wish,we can put some other work between
599: VecScatterBegin() and VecScatterEnd() to overlap the communication with
600: computation.
601: */
602: VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
603: VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
605: /*
606: Get pointers to vector data
607: */
608: VecGetArray(localX,&x);
609: VecGetArray(F,&f);
611: /*
612: Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
613: approximate one chosen for illustrative purpose only. Another point to notice
614: is that this is a local (completly parallel) calculation. In practical application
615: codes, function calculation time is a dominat portion of the overall execution time.
616: */
617: for (i=0; i < Nvlocal; i++) {
618: f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i];
619: for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]];
620: }
622: /*
623: Restore vectors
624: */
625: VecRestoreArray(localX,&x);
626: VecRestoreArray(F,&f);
627: /*VecView(F,PETSC_VIEWER_STDOUT_WORLD);*/
629: return 0;
630: }
632: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
633: /*
634: FormJacobian - Evaluates Jacobian matrix.
636: Input Parameters:
637: . snes - the SNES context
638: . X - input vector
639: . ptr - optional user-defined context, as set by SNESSetJacobian()
641: Output Parameters:
642: . A - Jacobian matrix
643: . B - optionally different preconditioning matrix
644: . flag - flag indicating matrix structure
646: */
647: PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr)
648: {
649: AppCtx *user = (AppCtx*)ptr;
650: PetscInt i,j,Nvlocal,col[50],ierr;
651: PetscScalar alpha,lambda,value[50];
652: Vec localX = user->localX;
653: VecScatter scatter;
654: PetscScalar *x;
655: #if defined(UNUSED_VARIABLES)
656: PetscScalar two = 2.0,one = 1.0;
657: PetscInt row,Nvglobal,Neglobal;
658: PetscInt *gloInd;
660: Nvglobal = user->Nvglobal;
661: Neglobal = user->Neglobal;
662: gloInd = user->gloInd;
663: #endif
665: /*printf("Entering into FormJacobian \n");*/
666: Nvlocal = user->Nvlocal;
667: lambda = user->non_lin_param;
668: alpha = user->lin_param;
669: scatter = user->scatter;
671: /*
672: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
673: described in the beginning of this code
675: First scatter the distributed vector X into local vector localX (that includes
676: values for ghost nodes. If we wish, we can put some other work between
677: VecScatterBegin() and VecScatterEnd() to overlap the communication with
678: computation.
679: */
680: VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
681: VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);
683: /*
684: Get pointer to vector data
685: */
686: VecGetArray(localX,&x);
688: for (i=0; i < Nvlocal; i++) {
689: col[0] = i;
690: value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha;
691: for (j = 0; j < user->itot[i]; j++) {
692: col[j+1] = user->AdjM[i][j];
693: value[j+1] = -1.0;
694: }
696: /*
697: Set the matrix values in the local ordering. Note that in order to use this
698: feature we must call the routine MatSetLocalToGlobalMapping() after the
699: matrix has been created.
700: */
701: MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);
702: }
704: /*
705: Assemble matrix, using the 2-step process:
706: MatAssemblyBegin(), MatAssemblyEnd().
707: Between these two calls, the pointer to vector data has been restored to
708: demonstrate the use of overlapping communicationn with computation.
709: */
710: MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);
711: VecRestoreArray(localX,&x);
712: MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);
714: /*
715: Tell the matrix we will never add a new nonzero location to the
716: matrix. If we do, it will generate an error.
717: */
718: MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
719: /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
720: return 0;
721: }