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