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>
12
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";
31
32 /* ------------------------------------------------------------------------
33
34 PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian.
35
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.
39
40 -----------------------------------------------------------------------*/
41
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 */
45
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;
63
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);
70
main(int argc,char ** argv)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;
104
105 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
106 Initialize program
107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108
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));
113
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.");
116
117 /*
118 Initialize problem parameters
119 */
120 user.Nvglobal = 16; /* Global # of vertices */
121 user.Neglobal = 18; /* Global # of elements */
122
123 PetscCall(PetscOptionsGetInt(NULL, NULL, "-vert", &user.Nvglobal, NULL));
124 PetscCall(PetscOptionsGetInt(NULL, NULL, "-elem", &user.Neglobal, NULL));
125
126 user.non_lin_param = 0.06;
127
128 PetscCall(PetscOptionsGetReal(NULL, NULL, "-nl_par", &user.non_lin_param, NULL));
129
130 user.lin_param = -1.0;
131
132 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lin_par", &user.lin_param, NULL));
133
134 user.Nvlocal = 0;
135 user.Nelocal = 0;
136
137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
138 Read the mesh and partitioning information
139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
140
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.
149
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");
158
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]));
174
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));
179
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)));
186
187 sscanf(str, form, &dtmp);
188 user.AdjM[user.Nvlocal][i] = dtmp;
189
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));
197
198 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
199 Create different orderings
200 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
201
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));
211
212 for (i = 0; i < user.Nvlocal; i++) pordering[i] = rstart + i;
213
214 /*
215 Create the AO object
216 */
217 PetscCall(AOCreateBasic(MPI_COMM_WORLD, user.Nvlocal, user.gloInd, pordering, &ao));
218 PetscCall(PetscFree(pordering));
219
220 /*
221 Keep the global indices for later use
222 */
223 PetscCall(PetscMalloc1(user.Nvlocal, &user.locInd));
224 PetscCall(PetscMalloc1(Nvneighborstotal, &tmp));
225
226 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
227 Demonstrate the use of AO functionality
228 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
229
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]));
233
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]));
242
243 tmp[j + jstart] = user.AdjM[i][j];
244 }
245 jstart += user.itot[i];
246 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
247 }
248
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));
255
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"));
259
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];
265
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 }
271
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
285
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;
294
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 }
302
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 }
315
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"));
320
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;
332
333 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]));
334 }
335 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
336 }
337
338 N = user.Nvglobal;
339
340 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
341 Create vector and matrix data structures
342 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
343
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));
353
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));
363
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));
373
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));
379
380 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
381 Create nonlinear solver context
382 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
383
384 PetscCall(SNESCreate(MPI_COMM_WORLD, &snes));
385 PetscCall(SNESSetType(snes, type));
386
387 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
388 Set routines for function and Jacobian evaluation
389 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
390 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user));
391
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;
398
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, (MatFDColoringFn *)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 }
415
416 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
417 Customize nonlinear solver; set runtime options
418 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
419
420 PetscCall(SNESSetFromOptions(snes));
421
422 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
423 Evaluate initial guess; then solve nonlinear system
424 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
425
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));
433
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));
440
441 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
442 Now solve the nonlinear system
443 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
444
445 PetscCall(SNESSolve(snes, NULL, x));
446 PetscCall(SNESGetIterationNumber(snes, &its));
447 PetscCall(SNESGetNonlinearStepFailures(snes, &nfails));
448
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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
454
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));
461
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 ----------------- */
485
486 /*
487 FormInitialGuess - Forms initial approximation.
488
489 Input Parameters:
490 user - user-defined application context
491 X - vector
492
493 Output Parameter:
494 X - vector
495 */
FormInitialGuess(AppCtx * user,Vec X)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;
505
506 Nvglobal = user->Nvglobal;
507 Neglobal = user->Neglobal;
508 lambda = user->non_lin_param;
509 alpha = user->lin_param;
510 #endif
511
512 PetscFunctionBeginUser;
513 Nvlocal = user->Nvlocal;
514 gloInd = user->gloInd;
515
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));
524
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];
529
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).
539
540 Input Parameters:
541 . snes - the SNES context
542 . X - input vector
543 . ptr - optional user-defined context, as set by SNESSetFunction()
544
545 Output Parameter:
546 . F - function vector
547 */
FormFunction(SNES snes,Vec X,Vec F,void * ptr)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;
562
563 Nvglobal = user->Nvglobal;
564 Neglobal = user->Neglobal;
565 gloInd = user->gloInd;
566 #endif
567
568 PetscFunctionBeginUser;
569 Nvlocal = user->Nvlocal;
570 lambda = user->non_lin_param;
571 alpha = user->lin_param;
572 scatter = user->scatter;
573
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
577
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));
585
586 /*
587 Get pointers to vector data
588 */
589 PetscCall(VecGetArray(localX, &x));
590 PetscCall(VecGetArray(F, &f));
591
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 }
602
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 }
611
612 /* -------------------- Evaluate Jacobian F'(x) -------------------- */
613 /*
614 FormJacobian - Evaluates Jacobian matrix.
615
616 Input Parameters:
617 . snes - the SNES context
618 . X - input vector
619 . ptr - optional user-defined context, as set by SNESSetJacobian()
620
621 Output Parameters:
622 . A - Jacobian matrix
623 . B - optionally different matrix used to construct the preconditioner
624
625 */
FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void * ptr)626 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
627 {
628 AppCtx *user = (AppCtx *)ptr;
629 PetscInt Nvlocal, col[50];
630 PetscScalar alpha, lambda, value[50];
631 Vec localX = user->localX;
632 VecScatter scatter;
633 PetscScalar *x;
634 #if defined(UNUSED_VARIABLES)
635 PetscScalar two = 2.0, one = 1.0;
636 PetscInt row, Nvglobal, Neglobal;
637 PetscInt *gloInd;
638
639 Nvglobal = user->Nvglobal;
640 Neglobal = user->Neglobal;
641 gloInd = user->gloInd;
642 #endif
643
644 PetscFunctionBeginUser;
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;
650
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
654
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 PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
661 PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
662
663 /*
664 Get pointer to vector data
665 */
666 PetscCall(VecGetArray(localX, &x));
667
668 for (PetscInt i = 0; i < Nvlocal; i++) {
669 col[0] = i;
670 value[0] = user->itot[i] - 2.0 * lambda * x[i] - alpha;
671 for (PetscInt j = 0; j < user->itot[i]; j++) {
672 col[j + 1] = user->AdjM[i][j];
673 value[j + 1] = -1.0;
674 }
675
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 PetscCall(MatSetValuesLocal(jac, 1, &i, 1 + user->itot[i], col, value, INSERT_VALUES));
682 }
683
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 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
691 PetscCall(VecRestoreArray(localX, &x));
692 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));
693
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 PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
699 /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
700 PetscFunctionReturn(PETSC_SUCCESS);
701 }
702
703 /*TEST
704
705 build:
706 requires: !complex
707
708 test:
709 nsize: 2
710 args: -snes_monitor_short
711 localrunfiles: options.inf adj.in
712
713 test:
714 suffix: 2
715 nsize: 2
716 args: -snes_monitor_short -fd_jacobian_coloring
717 localrunfiles: options.inf adj.in
718
719 TEST*/
720