Actual source code: ex10.c
1: /*
2: Include "petscsnes.h" so that we can use SNES solvers. Note that this
3: file automatically includes:
4: petscsys.h - base PETSc routines petscvec.h - vectors
5: petscmat.h - matrices
6: petscis.h - index sets petscksp.h - Krylov subspace methods
7: petscviewer.h - viewers petscpc.h - preconditioners
8: petscksp.h - linear solvers
9: */
10: #include <petscsnes.h>
11: #include <petscao.h>
13: static char help[] = "An Unstructured Grid Example.\n\
14: This example demonstrates how to solve a nonlinear system in parallel\n\
15: with SNES for an unstructured mesh. The mesh and partitioning information\n\
16: is read in an application defined ordering,which is later transformed\n\
17: into another convenient ordering (called the local ordering). The local\n\
18: ordering, apart from being efficient on cpu cycles and memory, allows\n\
19: the use of the SPMD model of parallel programming. After partitioning\n\
20: is done, scatters are created between local (sequential)and global\n\
21: (distributed) vectors. Finally, we set up the nonlinear solver context\n\
22: in the usual way as a structured grid (see\n\
23: petsc/src/snes/tutorials/ex5.c).\n\
24: This example also illustrates the use of parallel matrix coloring.\n\
25: The command line options include:\n\
26: -vert <Nv>, where Nv is the global number of nodes\n\
27: -elem <Ne>, where Ne is the global number of elements\n\
28: -nl_par <lambda>, where lambda is the multiplier for the non linear term (u*u) term\n\
29: -lin_par <alpha>, where alpha is the multiplier for the linear term (u)\n\
30: -fd_jacobian_coloring -mat_coloring_type lf\n";
32: /* ------------------------------------------------------------------------
34: PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian.
36: The Laplacian is approximated in the following way: each edge is given a weight
37: of one meaning that the diagonal term will have the weight equal to the degree
38: of a node. The off-diagonal terms will get a weight of -1.
40: -----------------------------------------------------------------------*/
42: #define MAX_ELEM 500 /* Maximum number of elements */
43: #define MAX_VERT 100 /* Maximum number of vertices */
44: #define MAX_VERT_ELEM 3 /* Vertices per element */
46: /*
47: Application-defined context for problem specific data
48: */
49: typedef struct {
50: PetscInt Nvglobal, Nvlocal; /* global and local number of vertices */
51: PetscInt Neglobal, Nelocal; /* global and local number of vertices */
52: PetscInt AdjM[MAX_VERT][50]; /* adjacency list of a vertex */
53: PetscInt itot[MAX_VERT]; /* total number of neighbors for a vertex */
54: PetscInt icv[MAX_ELEM][MAX_VERT_ELEM]; /* vertices belonging to an element */
55: PetscInt v2p[MAX_VERT]; /* processor number for a vertex */
56: PetscInt *locInd, *gloInd; /* local and global orderings for a node */
57: Vec localX, localF; /* local solution (u) and f(u) vectors */
58: PetscReal non_lin_param; /* nonlinear parameter for the PDE */
59: PetscReal lin_param; /* linear parameter for the PDE */
60: VecScatter scatter; /* scatter context for the local and
61: distributed vectors */
62: } AppCtx;
64: /*
65: User-defined routines
66: */
67: PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
68: PetscErrorCode FormFunction(SNES, Vec, Vec, void *);
69: PetscErrorCode FormInitialGuess(AppCtx *, Vec);
71: int main(int argc, char **argv)
72: {
73: SNES snes; /* SNES context */
74: SNESType type = SNESNEWTONLS; /* default nonlinear solution method */
75: Vec x, r; /* solution, residual vectors */
76: Mat Jac; /* Jacobian matrix */
77: AppCtx user; /* user-defined application context */
78: AO ao; /* Application Ordering object */
79: IS isglobal, islocal; /* global and local index sets */
80: PetscMPIInt rank, size; /* rank of a process, number of processors */
81: PetscInt rstart; /* starting index of PETSc ordering for a processor */
82: PetscInt nfails; /* number of unsuccessful Newton steps */
83: PetscInt bs = 1; /* block size for multicomponent systems */
84: PetscInt nvertices; /* number of local plus ghost nodes of a processor */
85: PetscInt *pordering; /* PETSc ordering */
86: PetscInt *vertices; /* list of all vertices (incl. ghost ones) on a processor */
87: PetscInt *verticesmask;
88: PetscInt *tmp;
89: PetscInt i, j, jstart, inode, nb, nbrs, Nvneighborstotal = 0;
90: PetscInt its, N;
91: PetscScalar *xx;
92: char str[256], form[256], part_name[256];
93: FILE *fptr, *fptr1;
94: ISLocalToGlobalMapping isl2g;
95: int dtmp;
96: #if defined(UNUSED_VARIABLES)
97: PetscDraw draw; /* drawing context */
98: PetscScalar *ff, *gg;
99: PetscReal tiny = 1.0e-10, zero = 0.0, one = 1.0, big = 1.0e+10;
100: PetscInt *tmp1, *tmp2;
101: #endif
102: MatFDColoring matfdcoloring = 0;
103: PetscBool fd_jacobian_coloring = PETSC_FALSE;
105: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106: Initialize program
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: PetscFunctionBeginUser;
110: PetscCall(PetscInitialize(&argc, &argv, "options.inf", help));
111: PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
112: PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));
114: /* The current input file options.inf is for 2 proc run only */
115: PetscCheck(size == 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This example currently runs on 2 procs only.");
117: /*
118: Initialize problem parameters
119: */
120: user.Nvglobal = 16; /* Global # of vertices */
121: user.Neglobal = 18; /* Global # of elements */
123: PetscCall(PetscOptionsGetInt(NULL, NULL, "-vert", &user.Nvglobal, NULL));
124: PetscCall(PetscOptionsGetInt(NULL, NULL, "-elem", &user.Neglobal, NULL));
126: user.non_lin_param = 0.06;
128: PetscCall(PetscOptionsGetReal(NULL, NULL, "-nl_par", &user.non_lin_param, NULL));
130: user.lin_param = -1.0;
132: PetscCall(PetscOptionsGetReal(NULL, NULL, "-lin_par", &user.lin_param, NULL));
134: user.Nvlocal = 0;
135: user.Nelocal = 0;
137: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138: Read the mesh and partitioning information
139: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
141: /*
142: Read the mesh and partitioning information from 'adj.in'.
143: The file format is as follows.
144: For each line the first entry is the processor rank where the
145: current node belongs. The second entry is the number of
146: neighbors of a node. The rest of the line is the adjacency
147: list of a node. Currently this file is set up to work on two
148: processors.
150: This is not a very good example of reading input. In the future,
151: we will put an example that shows the style that should be
152: used in a real application, where partitioning will be done
153: dynamically by calling partitioning routines (at present, we have
154: a ready interface to ParMeTiS).
155: */
156: fptr = fopen("adj.in", "r");
157: PetscCheck(fptr, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open adj.in");
159: /*
160: Each processor writes to the file output.<rank> where rank is the
161: processor's rank.
162: */
163: snprintf(part_name, PETSC_STATIC_ARRAY_LENGTH(part_name), "output.%d", rank);
164: fptr1 = fopen(part_name, "w");
165: PetscCheck(fptr1, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could no open output file");
166: PetscCall(PetscMalloc1(user.Nvglobal, &user.gloInd));
167: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Rank is %d\n", rank));
168: for (inode = 0; inode < user.Nvglobal; inode++) {
169: PetscCheck(fgets(str, 256, fptr), PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "fgets read failed");
170: sscanf(str, "%d", &dtmp);
171: user.v2p[inode] = dtmp;
172: if (user.v2p[inode] == rank) {
173: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Node %" PetscInt_FMT " belongs to processor %" PetscInt_FMT "\n", inode, user.v2p[inode]));
175: user.gloInd[user.Nvlocal] = inode;
176: sscanf(str, "%*d %d", &dtmp);
177: nbrs = dtmp;
178: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Number of neighbors for the vertex %" PetscInt_FMT " is %" PetscInt_FMT "\n", inode, nbrs));
180: user.itot[user.Nvlocal] = nbrs;
181: Nvneighborstotal += nbrs;
182: for (i = 0; i < user.itot[user.Nvlocal]; i++) {
183: form[0] = '\0';
184: for (j = 0; j < i + 2; j++) PetscCall(PetscStrlcat(form, "%*d ", sizeof(form)));
185: PetscCall(PetscStrlcat(form, "%d", sizeof(form)));
187: sscanf(str, form, &dtmp);
188: user.AdjM[user.Nvlocal][i] = dtmp;
190: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[user.Nvlocal][i]));
191: }
192: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
193: user.Nvlocal++;
194: }
195: }
196: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Total # of Local Vertices is %" PetscInt_FMT " \n", user.Nvlocal));
198: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199: Create different orderings
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
202: /*
203: Create the local ordering list for vertices. First a list using the PETSc global
204: ordering is created. Then we use the AO object to get the PETSc-to-application and
205: application-to-PETSc mappings. Each vertex also gets a local index (stored in the
206: locInd array).
207: */
208: PetscCallMPI(MPI_Scan(&user.Nvlocal, &rstart, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
209: rstart -= user.Nvlocal;
210: PetscCall(PetscMalloc1(user.Nvlocal, &pordering));
212: for (i = 0; i < user.Nvlocal; i++) pordering[i] = rstart + i;
214: /*
215: Create the AO object
216: */
217: PetscCall(AOCreateBasic(MPI_COMM_WORLD, user.Nvlocal, user.gloInd, pordering, &ao));
218: PetscCall(PetscFree(pordering));
220: /*
221: Keep the global indices for later use
222: */
223: PetscCall(PetscMalloc1(user.Nvlocal, &user.locInd));
224: PetscCall(PetscMalloc1(Nvneighborstotal, &tmp));
226: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227: Demonstrate the use of AO functionality
228: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
230: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Before AOApplicationToPetsc, local indices are : \n"));
231: for (i = 0; i < user.Nvlocal; i++) {
232: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, " %" PetscInt_FMT " ", user.gloInd[i]));
234: user.locInd[i] = user.gloInd[i];
235: }
236: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
237: jstart = 0;
238: for (i = 0; i < user.Nvlocal; i++) {
239: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Neghbors of local vertex %" PetscInt_FMT " are : ", user.gloInd[i]));
240: for (j = 0; j < user.itot[i]; j++) {
241: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]));
243: tmp[j + jstart] = user.AdjM[i][j];
244: }
245: jstart += user.itot[i];
246: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
247: }
249: /*
250: Now map the vlocal and neighbor lists to the PETSc ordering
251: */
252: PetscCall(AOApplicationToPetsc(ao, user.Nvlocal, user.locInd));
253: PetscCall(AOApplicationToPetsc(ao, Nvneighborstotal, tmp));
254: PetscCall(AODestroy(&ao));
256: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "After AOApplicationToPetsc, local indices are : \n"));
257: for (i = 0; i < user.Nvlocal; i++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, " %" PetscInt_FMT " ", user.locInd[i]));
258: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
260: jstart = 0;
261: for (i = 0; i < user.Nvlocal; i++) {
262: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Neghbors of local vertex %" PetscInt_FMT " are : ", user.locInd[i]));
263: for (j = 0; j < user.itot[i]; j++) {
264: user.AdjM[i][j] = tmp[j + jstart];
266: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]));
267: }
268: jstart += user.itot[i];
269: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
270: }
272: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
273: Extract the ghost vertex information for each processor
274: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
275: /*
276: Next, we need to generate a list of vertices required for this processor
277: and a local numbering scheme for all vertices required on this processor.
278: vertices - integer array of all vertices needed on this processor in PETSc
279: global numbering; this list consists of first the "locally owned"
280: vertices followed by the ghost vertices.
281: verticesmask - integer array that for each global vertex lists its local
282: vertex number (in vertices) + 1. If the global vertex is not
283: represented on this processor, then the corresponding
284: entry in verticesmask is zero
286: Note: vertices and verticesmask are both Nvglobal in length; this may
287: sound terribly non-scalable, but in fact is not so bad for a reasonable
288: number of processors. Importantly, it allows us to use NO SEARCHING
289: in setting up the data structures.
290: */
291: PetscCall(PetscMalloc1(user.Nvglobal, &vertices));
292: PetscCall(PetscCalloc1(user.Nvglobal, &verticesmask));
293: nvertices = 0;
295: /*
296: First load "owned vertices" into list
297: */
298: for (i = 0; i < user.Nvlocal; i++) {
299: vertices[nvertices++] = user.locInd[i];
300: verticesmask[user.locInd[i]] = nvertices;
301: }
303: /*
304: Now load ghost vertices into list
305: */
306: for (i = 0; i < user.Nvlocal; i++) {
307: for (j = 0; j < user.itot[i]; j++) {
308: nb = user.AdjM[i][j];
309: if (!verticesmask[nb]) {
310: vertices[nvertices++] = nb;
311: verticesmask[nb] = nvertices;
312: }
313: }
314: }
316: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
317: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "The array vertices is :\n"));
318: for (i = 0; i < nvertices; i++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", vertices[i]));
319: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
321: /*
322: Map the vertices listed in the neighbors to the local numbering from
323: the global ordering that they contained initially.
324: */
325: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
326: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "After mapping neighbors in the local contiguous ordering\n"));
327: for (i = 0; i < user.Nvlocal; i++) {
328: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Neghbors of local vertex %" PetscInt_FMT " are :\n", i));
329: for (j = 0; j < user.itot[i]; j++) {
330: nb = user.AdjM[i][j];
331: user.AdjM[i][j] = verticesmask[nb] - 1;
333: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]));
334: }
335: PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
336: }
338: N = user.Nvglobal;
340: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
341: Create vector and matrix data structures
342: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
344: /*
345: Create vector data structures
346: */
347: PetscCall(VecCreate(MPI_COMM_WORLD, &x));
348: PetscCall(VecSetSizes(x, user.Nvlocal, N));
349: PetscCall(VecSetFromOptions(x));
350: PetscCall(VecDuplicate(x, &r));
351: PetscCall(VecCreateSeq(MPI_COMM_SELF, bs * nvertices, &user.localX));
352: PetscCall(VecDuplicate(user.localX, &user.localF));
354: /*
355: Create the scatter between the global representation and the
356: local representation
357: */
358: PetscCall(ISCreateStride(MPI_COMM_SELF, bs * nvertices, 0, 1, &islocal));
359: PetscCall(ISCreateBlock(MPI_COMM_SELF, bs, nvertices, vertices, PETSC_COPY_VALUES, &isglobal));
360: PetscCall(VecScatterCreate(x, isglobal, user.localX, islocal, &user.scatter));
361: PetscCall(ISDestroy(&isglobal));
362: PetscCall(ISDestroy(&islocal));
364: /*
365: Create matrix data structure; Just to keep the example simple, we have not done any
366: preallocation of memory for the matrix. In real application code with big matrices,
367: preallocation should always be done to expedite the matrix creation.
368: */
369: PetscCall(MatCreate(MPI_COMM_WORLD, &Jac));
370: PetscCall(MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N));
371: PetscCall(MatSetFromOptions(Jac));
372: PetscCall(MatSetUp(Jac));
374: /*
375: The following routine allows us to set the matrix values in local ordering
376: */
377: PetscCall(ISLocalToGlobalMappingCreate(MPI_COMM_SELF, bs, nvertices, vertices, PETSC_COPY_VALUES, &isl2g));
378: PetscCall(MatSetLocalToGlobalMapping(Jac, isl2g, isl2g));
380: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
381: Create nonlinear solver context
382: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
384: PetscCall(SNESCreate(MPI_COMM_WORLD, &snes));
385: PetscCall(SNESSetType(snes, type));
387: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388: Set routines for function and Jacobian evaluation
389: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390: PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user));
392: PetscCall(PetscOptionsGetBool(NULL, NULL, "-fd_jacobian_coloring", &fd_jacobian_coloring, 0));
393: if (!fd_jacobian_coloring) {
394: PetscCall(SNESSetJacobian(snes, Jac, Jac, FormJacobian, (void *)&user));
395: } else { /* Use matfdcoloring */
396: ISColoring iscoloring;
397: MatColoring mc;
399: /* Get the data structure of Jac */
400: PetscCall(FormJacobian(snes, x, Jac, Jac, &user));
401: /* Create coloring context */
402: PetscCall(MatColoringCreate(Jac, &mc));
403: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
404: PetscCall(MatColoringSetFromOptions(mc));
405: PetscCall(MatColoringApply(mc, &iscoloring));
406: PetscCall(MatColoringDestroy(&mc));
407: PetscCall(MatFDColoringCreate(Jac, iscoloring, &matfdcoloring));
408: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormFunction, &user));
409: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
410: PetscCall(MatFDColoringSetUp(Jac, iscoloring, matfdcoloring));
411: /* PetscCall(MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD)); */
412: PetscCall(SNESSetJacobian(snes, Jac, Jac, SNESComputeJacobianDefaultColor, matfdcoloring));
413: PetscCall(ISColoringDestroy(&iscoloring));
414: }
416: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
417: Customize nonlinear solver; set runtime options
418: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
420: PetscCall(SNESSetFromOptions(snes));
422: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
423: Evaluate initial guess; then solve nonlinear system
424: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
426: /*
427: Note: The user should initialize the vector, x, with the initial guess
428: for the nonlinear solver prior to calling SNESSolve(). In particular,
429: to employ an initial guess of zero, the user should explicitly set
430: this vector to zero by calling VecSet().
431: */
432: PetscCall(FormInitialGuess(&user, x));
434: /*
435: Print the initial guess
436: */
437: PetscCall(VecGetArray(x, &xx));
438: for (inode = 0; inode < user.Nvlocal; inode++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Initial Solution at node %" PetscInt_FMT " is %f \n", inode, (double)PetscRealPart(xx[inode])));
439: PetscCall(VecRestoreArray(x, &xx));
441: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
442: Now solve the nonlinear system
443: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
445: PetscCall(SNESSolve(snes, NULL, x));
446: PetscCall(SNESGetIterationNumber(snes, &its));
447: PetscCall(SNESGetNonlinearStepFailures(snes, &nfails));
449: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
450: Print the output : solution vector and other information
451: Each processor writes to the file output.<rank> where rank is the
452: processor's rank.
453: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
455: PetscCall(VecGetArray(x, &xx));
456: for (inode = 0; inode < user.Nvlocal; inode++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Solution at node %" PetscInt_FMT " is %f \n", inode, (double)PetscRealPart(xx[inode])));
457: PetscCall(VecRestoreArray(x, &xx));
458: fclose(fptr1);
459: PetscCall(PetscPrintf(MPI_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT ", ", its));
460: PetscCall(PetscPrintf(MPI_COMM_WORLD, "number of unsuccessful steps = %" PetscInt_FMT "\n", nfails));
462: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
463: Free work space. All PETSc objects should be destroyed when they
464: are no longer needed.
465: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
466: PetscCall(PetscFree(user.gloInd));
467: PetscCall(PetscFree(user.locInd));
468: PetscCall(PetscFree(vertices));
469: PetscCall(PetscFree(verticesmask));
470: PetscCall(PetscFree(tmp));
471: PetscCall(VecScatterDestroy(&user.scatter));
472: PetscCall(ISLocalToGlobalMappingDestroy(&isl2g));
473: PetscCall(VecDestroy(&x));
474: PetscCall(VecDestroy(&r));
475: PetscCall(VecDestroy(&user.localX));
476: PetscCall(VecDestroy(&user.localF));
477: PetscCall(SNESDestroy(&snes));
478: PetscCall(MatDestroy(&Jac));
479: /* PetscCall(PetscDrawDestroy(draw));*/
480: if (fd_jacobian_coloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
481: PetscCall(PetscFinalize());
482: return 0;
483: }
484: /* -------------------- Form initial approximation ----------------- */
486: /*
487: FormInitialGuess - Forms initial approximation.
489: Input Parameters:
490: user - user-defined application context
491: X - vector
493: Output Parameter:
494: X - vector
495: */
496: PetscErrorCode FormInitialGuess(AppCtx *user, Vec X)
497: {
498: PetscInt Nvlocal;
499: PetscInt *gloInd;
500: PetscScalar *x;
501: #if defined(UNUSED_VARIABLES)
502: PetscReal temp1, temp, hx, hy, hxdhy, hydhx, sc;
503: PetscInt Neglobal, Nvglobal, j, row;
504: PetscReal alpha, lambda;
506: Nvglobal = user->Nvglobal;
507: Neglobal = user->Neglobal;
508: lambda = user->non_lin_param;
509: alpha = user->lin_param;
510: #endif
512: PetscFunctionBeginUser;
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: PetscCall(VecGetArray(X, &x));
525: /*
526: Compute initial guess over the locally owned part of the grid
527: */
528: for (PetscInt i = 0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i];
530: /*
531: Restore vector
532: */
533: PetscCall(VecRestoreArray(X, &x));
534: PetscFunctionReturn(PETSC_SUCCESS);
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 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: PetscFunctionBeginUser;
569: Nvlocal = user->Nvlocal;
570: lambda = user->non_lin_param;
571: alpha = user->lin_param;
572: scatter = user->scatter;
574: /*
575: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
576: described in the beginning of this code
578: First scatter the distributed vector X into local vector localX (that includes
579: values for ghost nodes. If we wish,we can put some other work between
580: VecScatterBegin() and VecScatterEnd() to overlap the communication with
581: computation.
582: */
583: PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
584: PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
586: /*
587: Get pointers to vector data
588: */
589: PetscCall(VecGetArray(localX, &x));
590: PetscCall(VecGetArray(F, &f));
592: /*
593: Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
594: approximate one chosen for illustrative purpose only. Another point to notice
595: is that this is a local (completely parallel) calculation. In practical application
596: codes, function calculation time is a dominat portion of the overall execution time.
597: */
598: for (PetscInt i = 0; i < Nvlocal; i++) {
599: f[i] = (user->itot[i] - alpha) * x[i] - lambda * x[i] * x[i];
600: for (PetscInt j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]];
601: }
603: /*
604: Restore vectors
605: */
606: PetscCall(VecRestoreArray(localX, &x));
607: PetscCall(VecRestoreArray(F, &f));
608: /* PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD)); */
609: PetscFunctionReturn(PETSC_SUCCESS);
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 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: PetscFunctionBeginUser;
646: /*printf("Entering into FormJacobian \n");*/
647: Nvlocal = user->Nvlocal;
648: lambda = user->non_lin_param;
649: alpha = user->lin_param;
650: scatter = user->scatter;
652: /*
653: PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
654: described in the beginning of this code
656: First scatter the distributed vector X into local vector localX (that includes
657: values for ghost nodes. If we wish, we can put some other work between
658: VecScatterBegin() and VecScatterEnd() to overlap the communication with
659: computation.
660: */
661: PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
662: PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
664: /*
665: Get pointer to vector data
666: */
667: PetscCall(VecGetArray(localX, &x));
669: for (PetscInt i = 0; i < Nvlocal; i++) {
670: col[0] = i;
671: value[0] = user->itot[i] - 2.0 * lambda * x[i] - alpha;
672: for (PetscInt j = 0; j < user->itot[i]; j++) {
673: col[j + 1] = user->AdjM[i][j];
674: value[j + 1] = -1.0;
675: }
677: /*
678: Set the matrix values in the local ordering. Note that in order to use this
679: feature we must call the routine MatSetLocalToGlobalMapping() after the
680: matrix has been created.
681: */
682: PetscCall(MatSetValuesLocal(jac, 1, &i, 1 + user->itot[i], col, value, INSERT_VALUES));
683: }
685: /*
686: Assemble matrix, using the 2-step process:
687: MatAssemblyBegin(), MatAssemblyEnd().
688: Between these two calls, the pointer to vector data has been restored to
689: demonstrate the use of overlapping communicationn with computation.
690: */
691: PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
692: PetscCall(VecRestoreArray(localX, &x));
693: PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
695: /*
696: Tell the matrix we will never add a new nonzero location to the
697: matrix. If we do, it will generate an error.
698: */
699: PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
700: /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
701: PetscFunctionReturn(PETSC_SUCCESS);
702: }
704: /*TEST
706: build:
707: requires: !complex
709: test:
710: nsize: 2
711: args: -snes_monitor_short
712: localrunfiles: options.inf adj.in
714: test:
715: suffix: 2
716: nsize: 2
717: args: -snes_monitor_short -fd_jacobian_coloring
718: localrunfiles: options.inf adj.in
720: TEST*/