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: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
111: PetscInitialize(&argc, &argv, "options.inf", help);
112: MPI_Comm_rank(MPI_COMM_WORLD, &rank);
113: MPI_Comm_size(MPI_COMM_WORLD, &size);
115: /* The current input file options.inf is for 2 proc run only */
118: /*
119: Initialize problem parameters
120: */
121: user.Nvglobal = 16; /* Global # of vertices */
122: user.Neglobal = 18; /* Global # of elements */
124: PetscOptionsGetInt(NULL, NULL, "-vert", &user.Nvglobal, NULL);
125: PetscOptionsGetInt(NULL, NULL, "-elem", &user.Neglobal, NULL);
127: user.non_lin_param = 0.06;
129: PetscOptionsGetReal(NULL, NULL, "-nl_par", &user.non_lin_param, NULL);
131: user.lin_param = -1.0;
133: PetscOptionsGetReal(NULL, NULL, "-lin_par", &user.lin_param, NULL);
135: user.Nvlocal = 0;
136: user.Nelocal = 0;
138: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
139: Read the mesh and partitioning information
140: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
142: /*
143: Read the mesh and partitioning information from 'adj.in'.
144: The file format is as follows.
145: For each line the first entry is the processor rank where the
146: current node belongs. The second entry is the number of
147: neighbors of a node. The rest of the line is the adjacency
148: list of a node. Currently this file is set up to work on two
149: processors.
151: This is not a very good example of reading input. In the future,
152: we will put an example that shows the style that should be
153: used in a real application, where partitioning will be done
154: dynamically by calling partitioning routines (at present, we have
155: a ready interface to ParMeTiS).
156: */
157: fptr = fopen("adj.in", "r");
160: /*
161: Each processor writes to the file output.<rank> where rank is the
162: processor's rank.
163: */
164: sprintf(part_name, "output.%d", rank);
165: fptr1 = fopen(part_name, "w");
167: PetscMalloc1(user.Nvglobal, &user.gloInd);
168: PetscFPrintf(PETSC_COMM_SELF, fptr1, "Rank is %d\n", rank);
169: for (inode = 0; inode < user.Nvglobal; inode++) {
171: sscanf(str, "%d", &dtmp);
172: user.v2p[inode] = dtmp;
173: if (user.v2p[inode] == rank) {
174: PetscFPrintf(PETSC_COMM_SELF, fptr1, "Node %" PetscInt_FMT " belongs to processor %" PetscInt_FMT "\n", inode, user.v2p[inode]);
176: user.gloInd[user.Nvlocal] = inode;
177: sscanf(str, "%*d %d", &dtmp);
178: nbrs = dtmp;
179: PetscFPrintf(PETSC_COMM_SELF, fptr1, "Number of neighbors for the vertex %" PetscInt_FMT " is %" PetscInt_FMT "\n", inode, nbrs);
181: user.itot[user.Nvlocal] = nbrs;
182: Nvneighborstotal += nbrs;
183: for (i = 0; i < user.itot[user.Nvlocal]; i++) {
184: form[0] = '\0';
185: for (j = 0; j < i + 2; j++) PetscStrlcat(form, "%*d ", sizeof(form));
186: PetscStrlcat(form, "%d", sizeof(form));
188: sscanf(str, form, &dtmp);
189: user.AdjM[user.Nvlocal][i] = dtmp;
191: PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", 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 %" PetscInt_FMT " \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, " %" PetscInt_FMT " ", 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 %" PetscInt_FMT " are : ", user.gloInd[i]);
241: for (j = 0; j < user.itot[i]; j++) {
242: PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", 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++) PetscFPrintf(PETSC_COMM_SELF, fptr1, " %" PetscInt_FMT " ", user.locInd[i]);
259: PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n");
261: jstart = 0;
262: for (i = 0; i < user.Nvlocal; i++) {
263: PetscFPrintf(PETSC_COMM_SELF, fptr1, "Neghbors of local vertex %" PetscInt_FMT " are : ", user.locInd[i]);
264: for (j = 0; j < user.itot[i]; j++) {
265: user.AdjM[i][j] = tmp[j + jstart];
267: PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]);
268: }
269: jstart += user.itot[i];
270: PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n");
271: }
273: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: Extract the ghost vertex information for each processor
275: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
276: /*
277: Next, we need to generate a list of vertices required for this processor
278: and a local numbering scheme for all vertices required on this processor.
279: vertices - integer array of all vertices needed on this processor in PETSc
280: global numbering; this list consists of first the "locally owned"
281: vertices followed by the ghost vertices.
282: verticesmask - integer array that for each global vertex lists its local
283: vertex number (in vertices) + 1. If the global vertex is not
284: represented on this processor, then the corresponding
285: entry in verticesmask is zero
287: Note: vertices and verticesmask are both Nvglobal in length; this may
288: sound terribly non-scalable, but in fact is not so bad for a reasonable
289: number of processors. Importantly, it allows us to use NO SEARCHING
290: in setting up the data structures.
291: */
292: PetscMalloc1(user.Nvglobal, &vertices);
293: PetscCalloc1(user.Nvglobal, &verticesmask);
294: nvertices = 0;
296: /*
297: First load "owned vertices" into list
298: */
299: for (i = 0; i < user.Nvlocal; i++) {
300: vertices[nvertices++] = user.locInd[i];
301: verticesmask[user.locInd[i]] = nvertices;
302: }
304: /*
305: Now load ghost vertices into list
306: */
307: for (i = 0; i < user.Nvlocal; i++) {
308: for (j = 0; j < user.itot[i]; j++) {
309: nb = user.AdjM[i][j];
310: if (!verticesmask[nb]) {
311: vertices[nvertices++] = nb;
312: verticesmask[nb] = nvertices;
313: }
314: }
315: }
317: PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n");
318: PetscFPrintf(PETSC_COMM_SELF, fptr1, "The array vertices is :\n");
319: for (i = 0; i < nvertices; i++) PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", vertices[i]);
320: PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n");
322: /*
323: Map the vertices listed in the neighbors to the local numbering from
324: the global ordering that they contained initially.
325: */
326: PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n");
327: PetscFPrintf(PETSC_COMM_SELF, fptr1, "After mapping neighbors in the local contiguous ordering\n");
328: for (i = 0; i < user.Nvlocal; i++) {
329: PetscFPrintf(PETSC_COMM_SELF, fptr1, "Neghbors of local vertex %" PetscInt_FMT " are :\n", i);
330: for (j = 0; j < user.itot[i]; j++) {
331: nb = user.AdjM[i][j];
332: user.AdjM[i][j] = verticesmask[nb] - 1;
334: PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]);
335: }
336: PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n");
337: }
339: N = user.Nvglobal;
341: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
342: Create vector and matrix data structures
343: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
345: /*
346: Create vector data structures
347: */
348: VecCreate(MPI_COMM_WORLD, &x);
349: VecSetSizes(x, user.Nvlocal, N);
350: VecSetFromOptions(x);
351: VecDuplicate(x, &r);
352: VecCreateSeq(MPI_COMM_SELF, bs * nvertices, &user.localX);
353: VecDuplicate(user.localX, &user.localF);
355: /*
356: Create the scatter between the global representation and the
357: local representation
358: */
359: ISCreateStride(MPI_COMM_SELF, bs * nvertices, 0, 1, &islocal);
360: ISCreateBlock(MPI_COMM_SELF, bs, nvertices, vertices, PETSC_COPY_VALUES, &isglobal);
361: VecScatterCreate(x, isglobal, user.localX, islocal, &user.scatter);
362: ISDestroy(&isglobal);
363: ISDestroy(&islocal);
365: /*
366: Create matrix data structure; Just to keep the example simple, we have not done any
367: preallocation of memory for the matrix. In real application code with big matrices,
368: preallocation should always be done to expedite the matrix creation.
369: */
370: MatCreate(MPI_COMM_WORLD, &Jac);
371: MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N);
372: MatSetFromOptions(Jac);
373: MatSetUp(Jac);
375: /*
376: The following routine allows us to set the matrix values in local ordering
377: */
378: ISLocalToGlobalMappingCreate(MPI_COMM_SELF, bs, nvertices, vertices, PETSC_COPY_VALUES, &isl2g);
379: MatSetLocalToGlobalMapping(Jac, isl2g, isl2g);
381: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
382: Create nonlinear solver context
383: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
385: SNESCreate(MPI_COMM_WORLD, &snes);
386: SNESSetType(snes, type);
388: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
389: Set routines for function and Jacobian evaluation
390: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391: SNESSetFunction(snes, r, FormFunction, (void *)&user);
393: PetscOptionsGetBool(NULL, NULL, "-fd_jacobian_coloring", &fd_jacobian_coloring, 0);
394: if (!fd_jacobian_coloring) {
395: SNESSetJacobian(snes, Jac, Jac, FormJacobian, (void *)&user);
396: } else { /* Use matfdcoloring */
397: ISColoring iscoloring;
398: MatColoring mc;
400: /* Get the data structure of Jac */
401: FormJacobian(snes, x, Jac, Jac, &user);
402: /* Create coloring context */
403: MatColoringCreate(Jac, &mc);
404: MatColoringSetType(mc, MATCOLORINGSL);
405: MatColoringSetFromOptions(mc);
406: MatColoringApply(mc, &iscoloring);
407: MatColoringDestroy(&mc);
408: MatFDColoringCreate(Jac, iscoloring, &matfdcoloring);
409: MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormFunction, &user);
410: MatFDColoringSetFromOptions(matfdcoloring);
411: MatFDColoringSetUp(Jac, iscoloring, matfdcoloring);
412: /* MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD); */
413: SNESSetJacobian(snes, Jac, Jac, SNESComputeJacobianDefaultColor, matfdcoloring);
414: ISColoringDestroy(&iscoloring);
415: }
417: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
418: Customize nonlinear solver; set runtime options
419: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
421: SNESSetFromOptions(snes);
423: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
424: Evaluate initial guess; then solve nonlinear system
425: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
427: /*
428: Note: The user should initialize the vector, x, with the initial guess
429: for the nonlinear solver prior to calling SNESSolve(). In particular,
430: to employ an initial guess of zero, the user should explicitly set
431: this vector to zero by calling VecSet().
432: */
433: FormInitialGuess(&user, x);
435: /*
436: Print the initial guess
437: */
438: VecGetArray(x, &xx);
439: for (inode = 0; inode < user.Nvlocal; inode++) PetscFPrintf(PETSC_COMM_SELF, fptr1, "Initial Solution at node %" PetscInt_FMT " is %f \n", inode, (double)PetscRealPart(xx[inode]));
440: VecRestoreArray(x, &xx);
442: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
443: Now solve the nonlinear system
444: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
446: SNESSolve(snes, NULL, x);
447: SNESGetIterationNumber(snes, &its);
448: SNESGetNonlinearStepFailures(snes, &nfails);
450: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
451: Print the output : solution vector and other information
452: Each processor writes to the file output.<rank> where rank is the
453: processor's rank.
454: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
456: VecGetArray(x, &xx);
457: for (inode = 0; inode < user.Nvlocal; inode++) PetscFPrintf(PETSC_COMM_SELF, fptr1, "Solution at node %" PetscInt_FMT " is %f \n", inode, (double)PetscRealPart(xx[inode]));
458: VecRestoreArray(x, &xx);
459: fclose(fptr1);
460: PetscPrintf(MPI_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT ", ", its);
461: PetscPrintf(MPI_COMM_WORLD, "number of unsuccessful steps = %" PetscInt_FMT "\n", nfails);
463: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
464: Free work space. All PETSc objects should be destroyed when they
465: are no longer needed.
466: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
467: PetscFree(user.gloInd);
468: PetscFree(user.locInd);
469: PetscFree(vertices);
470: PetscFree(verticesmask);
471: PetscFree(tmp);
472: VecScatterDestroy(&user.scatter);
473: ISLocalToGlobalMappingDestroy(&isl2g);
474: VecDestroy(&x);
475: VecDestroy(&r);
476: VecDestroy(&user.localX);
477: VecDestroy(&user.localF);
478: SNESDestroy(&snes);
479: MatDestroy(&Jac);
480: /* PetscDrawDestroy(draw);*/
481: if (fd_jacobian_coloring) MatFDColoringDestroy(&matfdcoloring);
482: PetscFinalize();
483: return 0;
484: }
485: /* -------------------- Form initial approximation ----------------- */
487: /*
488: FormInitialGuess - Forms initial approximation.
490: Input Parameters:
491: user - user-defined application context
492: X - vector
494: Output Parameter:
495: X - vector
496: */
497: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
498: {
499: PetscInt i, Nvlocal;
500: PetscInt *gloInd;
501: PetscScalar *x;
502: #if defined(UNUSED_VARIABLES)
503: PetscReal temp1, temp, hx, hy, hxdhy, hydhx, sc;
504: PetscInt Neglobal, Nvglobal, j, row;
505: PetscReal alpha, lambda;
507: Nvglobal = user->Nvglobal;
508: Neglobal = user->Neglobal;
509: lambda = user->non_lin_param;
510: alpha = user->lin_param;
511: #endif
513: Nvlocal = user->Nvlocal;
514: gloInd = user->gloInd;
516: /*
517: Get a pointer to vector data.
518: - For default PETSc vectors, VecGetArray() returns a pointer to
519: the data array. Otherwise, the routine is implementation dependent.
520: - You MUST call VecRestoreArray() when you no longer need access to
521: the array.
522: */
523: VecGetArray(X, &x);
525: /*
526: Compute initial guess over the locally owned part of the grid
527: */
528: for (i = 0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i];
530: /*
531: Restore vector
532: */
533: VecRestoreArray(X, &x);
534: return 0;
535: }
536: /* -------------------- Evaluate Function F(x) --------------------- */
537: /*
538: FormFunction - Evaluates nonlinear function, F(x).
540: Input Parameters:
541: . snes - the SNES context
542: . X - input vector
543: . ptr - optional user-defined context, as set by SNESSetFunction()
545: Output Parameter:
546: . F - function vector
547: */
548: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
549: {
550: AppCtx *user = (AppCtx *)ptr;
551: PetscInt i, j, Nvlocal;
552: PetscReal alpha, lambda;
553: PetscScalar *x, *f;
554: VecScatter scatter;
555: Vec localX = user->localX;
556: #if defined(UNUSED_VARIABLES)
557: PetscScalar ut, ub, ul, ur, u, *g, sc, uyy, uxx;
558: PetscReal hx, hy, hxdhy, hydhx;
559: PetscReal two = 2.0, one = 1.0;
560: PetscInt Nvglobal, Neglobal, row;
561: PetscInt *gloInd;
563: Nvglobal = user->Nvglobal;
564: Neglobal = user->Neglobal;
565: gloInd = user->gloInd;
566: #endif
568: Nvlocal = user->Nvlocal;
569: lambda = user->non_lin_param;
570: alpha = user->lin_param;
571: scatter = user->scatter;
573: /*
574: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
575: described in the beginning of this code
577: First scatter the distributed vector X into local vector localX (that includes
578: values for ghost nodes. If we wish,we can put some other work between
579: VecScatterBegin() and VecScatterEnd() to overlap the communication with
580: computation.
581: */
582: VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD);
583: VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD);
585: /*
586: Get pointers to vector data
587: */
588: VecGetArray(localX, &x);
589: VecGetArray(F, &f);
591: /*
592: Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
593: approximate one chosen for illustrative purpose only. Another point to notice
594: is that this is a local (completely parallel) calculation. In practical application
595: codes, function calculation time is a dominat portion of the overall execution time.
596: */
597: for (i = 0; i < Nvlocal; i++) {
598: f[i] = (user->itot[i] - alpha) * x[i] - lambda * x[i] * x[i];
599: for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]];
600: }
602: /*
603: Restore vectors
604: */
605: VecRestoreArray(localX, &x);
606: VecRestoreArray(F, &f);
607: /* VecView(F,PETSC_VIEWER_STDOUT_WORLD); */
609: return 0;
610: }
612: /* -------------------- Evaluate Jacobian F'(x) -------------------- */
613: /*
614: FormJacobian - Evaluates Jacobian matrix.
616: Input Parameters:
617: . snes - the SNES context
618: . X - input vector
619: . ptr - optional user-defined context, as set by SNESSetJacobian()
621: Output Parameters:
622: . A - Jacobian matrix
623: . B - optionally different preconditioning matrix
624: . flag - flag indicating matrix structure
626: */
627: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
628: {
629: AppCtx *user = (AppCtx *)ptr;
630: PetscInt i, j, Nvlocal, col[50];
631: PetscScalar alpha, lambda, value[50];
632: Vec localX = user->localX;
633: VecScatter scatter;
634: PetscScalar *x;
635: #if defined(UNUSED_VARIABLES)
636: PetscScalar two = 2.0, one = 1.0;
637: PetscInt row, Nvglobal, Neglobal;
638: PetscInt *gloInd;
640: Nvglobal = user->Nvglobal;
641: Neglobal = user->Neglobal;
642: gloInd = user->gloInd;
643: #endif
645: /*printf("Entering into FormJacobian \n");*/
646: Nvlocal = user->Nvlocal;
647: lambda = user->non_lin_param;
648: alpha = user->lin_param;
649: scatter = user->scatter;
651: /*
652: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
653: described in the beginning of this code
655: First scatter the distributed vector X into local vector localX (that includes
656: values for ghost nodes. If we wish, we can put some other work between
657: VecScatterBegin() and VecScatterEnd() to overlap the communication with
658: computation.
659: */
660: VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD);
661: VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD);
663: /*
664: Get pointer to vector data
665: */
666: VecGetArray(localX, &x);
668: for (i = 0; i < Nvlocal; i++) {
669: col[0] = i;
670: value[0] = user->itot[i] - 2.0 * lambda * x[i] - alpha;
671: for (j = 0; j < user->itot[i]; j++) {
672: col[j + 1] = user->AdjM[i][j];
673: value[j + 1] = -1.0;
674: }
676: /*
677: Set the matrix values in the local ordering. Note that in order to use this
678: feature we must call the routine MatSetLocalToGlobalMapping() after the
679: matrix has been created.
680: */
681: MatSetValuesLocal(jac, 1, &i, 1 + user->itot[i], col, value, INSERT_VALUES);
682: }
684: /*
685: Assemble matrix, using the 2-step process:
686: MatAssemblyBegin(), MatAssemblyEnd().
687: Between these two calls, the pointer to vector data has been restored to
688: demonstrate the use of overlapping communicationn with computation.
689: */
690: MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY);
691: VecRestoreArray(localX, &x);
692: MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY);
694: /*
695: Tell the matrix we will never add a new nonzero location to the
696: matrix. If we do, it will generate an error.
697: */
698: MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE);
699: /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
700: return 0;
701: }
703: /*TEST
705: build:
706: requires: !complex
708: test:
709: nsize: 2
710: args: -snes_monitor_short
711: localrunfiles: options.inf adj.in
713: test:
714: suffix: 2
715: nsize: 2
716: args: -snes_monitor_short -fd_jacobian_coloring
717: localrunfiles: options.inf adj.in
719: TEST*/