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