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:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

110:   PetscFunctionBeginUser;
111:   PetscCall(PetscInitialize(&argc, &argv, "options.inf", help));
112:   PetscCallMPI(MPI_Comm_rank(MPI_COMM_WORLD, &rank));
113:   PetscCallMPI(MPI_Comm_size(MPI_COMM_WORLD, &size));

115:   /* The current input file options.inf is for 2 proc run only */
116:   PetscCheck(size == 2, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "This example currently runs on 2 procs only.");

118:   /*
119:      Initialize problem parameters
120:   */
121:   user.Nvglobal = 16; /* Global # of vertices  */
122:   user.Neglobal = 18; /* Global # of elements  */

124:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-vert", &user.Nvglobal, NULL));
125:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-elem", &user.Neglobal, NULL));

127:   user.non_lin_param = 0.06;

129:   PetscCall(PetscOptionsGetReal(NULL, NULL, "-nl_par", &user.non_lin_param, NULL));

131:   user.lin_param = -1.0;

133:   PetscCall(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");
158:   PetscCheck(fptr, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could not open adj.in");

160:   /*
161:      Each processor writes to the file output.<rank> where rank is the
162:      processor's rank.
163:   */
164:   snprintf(part_name, sizeof(part_name) / sizeof(part_name[0]), "output.%d", rank);
165:   fptr1 = fopen(part_name, "w");
166:   PetscCheck(fptr1, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Could no open output file");
167:   PetscCall(PetscMalloc1(user.Nvglobal, &user.gloInd));
168:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Rank is %d\n", rank));
169:   for (inode = 0; inode < user.Nvglobal; inode++) {
170:     PetscCheck(fgets(str, 256, fptr), PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "fgets read failed");
171:     sscanf(str, "%d", &dtmp);
172:     user.v2p[inode] = dtmp;
173:     if (user.v2p[inode] == rank) {
174:       PetscCall(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:       PetscCall(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++) PetscCall(PetscStrlcat(form, "%*d ", sizeof(form)));
186:         PetscCall(PetscStrlcat(form, "%d", sizeof(form)));

188:         sscanf(str, form, &dtmp);
189:         user.AdjM[user.Nvlocal][i] = dtmp;

191:         PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[user.Nvlocal][i]));
192:       }
193:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
194:       user.Nvlocal++;
195:     }
196:   }
197:   PetscCall(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:   PetscCallMPI(MPI_Scan(&user.Nvlocal, &rstart, 1, MPIU_INT, MPI_SUM, PETSC_COMM_WORLD));
210:   rstart -= user.Nvlocal;
211:   PetscCall(PetscMalloc1(user.Nvlocal, &pordering));

213:   for (i = 0; i < user.Nvlocal; i++) pordering[i] = rstart + i;

215:   /*
216:     Create the AO object
217:   */
218:   PetscCall(AOCreateBasic(MPI_COMM_WORLD, user.Nvlocal, user.gloInd, pordering, &ao));
219:   PetscCall(PetscFree(pordering));

221:   /*
222:     Keep the global indices for later use
223:   */
224:   PetscCall(PetscMalloc1(user.Nvlocal, &user.locInd));
225:   PetscCall(PetscMalloc1(Nvneighborstotal, &tmp));

227:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228:     Demonstrate the use of AO functionality
229:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

231:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "Before AOApplicationToPetsc, local indices are : \n"));
232:   for (i = 0; i < user.Nvlocal; i++) {
233:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, " %" PetscInt_FMT " ", user.gloInd[i]));

235:     user.locInd[i] = user.gloInd[i];
236:   }
237:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
238:   jstart = 0;
239:   for (i = 0; i < user.Nvlocal; i++) {
240:     PetscCall(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:       PetscCall(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:     PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
248:   }

250:   /*
251:     Now map the vlocal and neighbor lists to the PETSc ordering
252:   */
253:   PetscCall(AOApplicationToPetsc(ao, user.Nvlocal, user.locInd));
254:   PetscCall(AOApplicationToPetsc(ao, Nvneighborstotal, tmp));
255:   PetscCall(AODestroy(&ao));

257:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "After AOApplicationToPetsc, local indices are : \n"));
258:   for (i = 0; i < user.Nvlocal; i++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, " %" PetscInt_FMT " ", user.locInd[i]));
259:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));

261:   jstart = 0;
262:   for (i = 0; i < user.Nvlocal; i++) {
263:     PetscCall(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:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]));
268:     }
269:     jstart += user.itot[i];
270:     PetscCall(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:   PetscCall(PetscMalloc1(user.Nvglobal, &vertices));
293:   PetscCall(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:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
318:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "The array vertices is :\n"));
319:   for (i = 0; i < nvertices; i++) PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", vertices[i]));
320:   PetscCall(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:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n"));
327:   PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "After mapping neighbors in the local contiguous ordering\n"));
328:   for (i = 0; i < user.Nvlocal; i++) {
329:     PetscCall(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:       PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j]));
335:     }
336:     PetscCall(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:   PetscCall(VecCreate(MPI_COMM_WORLD, &x));
349:   PetscCall(VecSetSizes(x, user.Nvlocal, N));
350:   PetscCall(VecSetFromOptions(x));
351:   PetscCall(VecDuplicate(x, &r));
352:   PetscCall(VecCreateSeq(MPI_COMM_SELF, bs * nvertices, &user.localX));
353:   PetscCall(VecDuplicate(user.localX, &user.localF));

355:   /*
356:     Create the scatter between the global representation and the
357:     local representation
358:   */
359:   PetscCall(ISCreateStride(MPI_COMM_SELF, bs * nvertices, 0, 1, &islocal));
360:   PetscCall(ISCreateBlock(MPI_COMM_SELF, bs, nvertices, vertices, PETSC_COPY_VALUES, &isglobal));
361:   PetscCall(VecScatterCreate(x, isglobal, user.localX, islocal, &user.scatter));
362:   PetscCall(ISDestroy(&isglobal));
363:   PetscCall(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:   PetscCall(MatCreate(MPI_COMM_WORLD, &Jac));
371:   PetscCall(MatSetSizes(Jac, PETSC_DECIDE, PETSC_DECIDE, N, N));
372:   PetscCall(MatSetFromOptions(Jac));
373:   PetscCall(MatSetUp(Jac));

375:   /*
376:     The following routine allows us to set the matrix values in local ordering
377:   */
378:   PetscCall(ISLocalToGlobalMappingCreate(MPI_COMM_SELF, bs, nvertices, vertices, PETSC_COPY_VALUES, &isl2g));
379:   PetscCall(MatSetLocalToGlobalMapping(Jac, isl2g, isl2g));

381:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
382:      Create nonlinear solver context
383:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

385:   PetscCall(SNESCreate(MPI_COMM_WORLD, &snes));
386:   PetscCall(SNESSetType(snes, type));

388:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
389:     Set routines for function and Jacobian evaluation
390:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
391:   PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user));

393:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-fd_jacobian_coloring", &fd_jacobian_coloring, 0));
394:   if (!fd_jacobian_coloring) {
395:     PetscCall(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:     PetscCall(FormJacobian(snes, x, Jac, Jac, &user));
402:     /* Create coloring context */
403:     PetscCall(MatColoringCreate(Jac, &mc));
404:     PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
405:     PetscCall(MatColoringSetFromOptions(mc));
406:     PetscCall(MatColoringApply(mc, &iscoloring));
407:     PetscCall(MatColoringDestroy(&mc));
408:     PetscCall(MatFDColoringCreate(Jac, iscoloring, &matfdcoloring));
409:     PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))FormFunction, &user));
410:     PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
411:     PetscCall(MatFDColoringSetUp(Jac, iscoloring, matfdcoloring));
412:     /* PetscCall(MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD)); */
413:     PetscCall(SNESSetJacobian(snes, Jac, Jac, SNESComputeJacobianDefaultColor, matfdcoloring));
414:     PetscCall(ISColoringDestroy(&iscoloring));
415:   }

417:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
418:      Customize nonlinear solver; set runtime options
419:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

421:   PetscCall(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:   PetscCall(FormInitialGuess(&user, x));

435:   /*
436:     Print the initial guess
437:   */
438:   PetscCall(VecGetArray(x, &xx));
439:   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])));
440:   PetscCall(VecRestoreArray(x, &xx));

442:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
443:      Now solve the nonlinear system
444:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

446:   PetscCall(SNESSolve(snes, NULL, x));
447:   PetscCall(SNESGetIterationNumber(snes, &its));
448:   PetscCall(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:   PetscCall(VecGetArray(x, &xx));
457:   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])));
458:   PetscCall(VecRestoreArray(x, &xx));
459:   fclose(fptr1);
460:   PetscCall(PetscPrintf(MPI_COMM_WORLD, "number of SNES iterations = %" PetscInt_FMT ", ", its));
461:   PetscCall(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:   PetscCall(PetscFree(user.gloInd));
468:   PetscCall(PetscFree(user.locInd));
469:   PetscCall(PetscFree(vertices));
470:   PetscCall(PetscFree(verticesmask));
471:   PetscCall(PetscFree(tmp));
472:   PetscCall(VecScatterDestroy(&user.scatter));
473:   PetscCall(ISLocalToGlobalMappingDestroy(&isl2g));
474:   PetscCall(VecDestroy(&x));
475:   PetscCall(VecDestroy(&r));
476:   PetscCall(VecDestroy(&user.localX));
477:   PetscCall(VecDestroy(&user.localF));
478:   PetscCall(SNESDestroy(&snes));
479:   PetscCall(MatDestroy(&Jac));
480:   /* PetscCall(PetscDrawDestroy(draw));*/
481:   if (fd_jacobian_coloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
482:   PetscCall(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     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:   PetscFunctionBeginUser;
514:   Nvlocal = user->Nvlocal;
515:   gloInd  = user->gloInd;

517:   /*
518:      Get a pointer to vector data.
519:        - For default PETSc vectors, VecGetArray() returns a pointer to
520:          the data array.  Otherwise, the routine is implementation dependent.
521:        - You MUST call VecRestoreArray() when you no longer need access to
522:          the array.
523:   */
524:   PetscCall(VecGetArray(X, &x));

526:   /*
527:      Compute initial guess over the locally owned part of the grid
528:   */
529:   for (PetscInt i = 0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i];

531:   /*
532:      Restore vector
533:   */
534:   PetscCall(VecRestoreArray(X, &x));
535:   PetscFunctionReturn(PETSC_SUCCESS);
536: }
537: /* --------------------  Evaluate Function F(x) --------------------- */
538: /*
539:    FormFunction - Evaluates nonlinear function, F(x).

541:    Input Parameters:
542: .  snes - the SNES context
543: .  X - input vector
544: .  ptr - optional user-defined context, as set by SNESSetFunction()

546:    Output Parameter:
547: .  F - function vector
548:  */
549: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr)
550: {
551:   AppCtx      *user = (AppCtx *)ptr;
552:   PetscInt     Nvlocal;
553:   PetscReal    alpha, lambda;
554:   PetscScalar *x, *f;
555:   VecScatter   scatter;
556:   Vec          localX = user->localX;
557: #if defined(UNUSED_VARIABLES)
558:   PetscScalar ut, ub, ul, ur, u, *g, sc, uyy, uxx;
559:   PetscReal   hx, hy, hxdhy, hydhx;
560:   PetscReal   two = 2.0, one = 1.0;
561:   PetscInt    Nvglobal, Neglobal, row;
562:   PetscInt   *gloInd;

564:   Nvglobal = user->Nvglobal;
565:   Neglobal = user->Neglobal;
566:   gloInd   = user->gloInd;
567: #endif

569:   PetscFunctionBeginUser;
570:   Nvlocal = user->Nvlocal;
571:   lambda  = user->non_lin_param;
572:   alpha   = user->lin_param;
573:   scatter = user->scatter;

575:   /*
576:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
577:      described in the beginning of this code

579:      First scatter the distributed vector X into local vector localX (that includes
580:      values for ghost nodes. If we wish,we can put some other work between
581:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
582:      computation.
583:  */
584:   PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
585:   PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));

587:   /*
588:      Get pointers to vector data
589:   */
590:   PetscCall(VecGetArray(localX, &x));
591:   PetscCall(VecGetArray(F, &f));

593:   /*
594:     Now compute the f(x). As mentioned earlier, the computed Laplacian is just an
595:     approximate one chosen for illustrative purpose only. Another point to notice
596:     is that this is a local (completely parallel) calculation. In practical application
597:     codes, function calculation time is a dominat portion of the overall execution time.
598:   */
599:   for (PetscInt i = 0; i < Nvlocal; i++) {
600:     f[i] = (user->itot[i] - alpha) * x[i] - lambda * x[i] * x[i];
601:     for (PetscInt j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]];
602:   }

604:   /*
605:      Restore vectors
606:   */
607:   PetscCall(VecRestoreArray(localX, &x));
608:   PetscCall(VecRestoreArray(F, &f));
609:   /* PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD)); */

611:   PetscFunctionReturn(PETSC_SUCCESS);
612: }

614: /* --------------------  Evaluate Jacobian F'(x) -------------------- */
615: /*
616:    FormJacobian - Evaluates Jacobian matrix.

618:    Input Parameters:
619: .  snes - the SNES context
620: .  X - input vector
621: .  ptr - optional user-defined context, as set by SNESSetJacobian()

623:    Output Parameters:
624: .  A - Jacobian matrix
625: .  B - optionally different preconditioning matrix
626: .  flag - flag indicating matrix structure

628: */
629: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr)
630: {
631:   AppCtx      *user = (AppCtx *)ptr;
632:   PetscInt     Nvlocal, col[50];
633:   PetscScalar  alpha, lambda, value[50];
634:   Vec          localX = user->localX;
635:   VecScatter   scatter;
636:   PetscScalar *x;
637: #if defined(UNUSED_VARIABLES)
638:   PetscScalar two = 2.0, one = 1.0;
639:   PetscInt    row, Nvglobal, Neglobal;
640:   PetscInt   *gloInd;

642:   Nvglobal = user->Nvglobal;
643:   Neglobal = user->Neglobal;
644:   gloInd   = user->gloInd;
645: #endif

647:   PetscFunctionBeginUser;
648:   /*printf("Entering into FormJacobian \n");*/
649:   Nvlocal = user->Nvlocal;
650:   lambda  = user->non_lin_param;
651:   alpha   = user->lin_param;
652:   scatter = user->scatter;

654:   /*
655:      PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as
656:      described in the beginning of this code

658:      First scatter the distributed vector X into local vector localX (that includes
659:      values for ghost nodes. If we wish, we can put some other work between
660:      VecScatterBegin() and VecScatterEnd() to overlap the communication with
661:      computation.
662:   */
663:   PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
664:   PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));

666:   /*
667:      Get pointer to vector data
668:   */
669:   PetscCall(VecGetArray(localX, &x));

671:   for (PetscInt i = 0; i < Nvlocal; i++) {
672:     col[0]   = i;
673:     value[0] = user->itot[i] - 2.0 * lambda * x[i] - alpha;
674:     for (PetscInt j = 0; j < user->itot[i]; j++) {
675:       col[j + 1]   = user->AdjM[i][j];
676:       value[j + 1] = -1.0;
677:     }

679:     /*
680:       Set the matrix values in the local ordering. Note that in order to use this
681:       feature we must call the routine MatSetLocalToGlobalMapping() after the
682:       matrix has been created.
683:     */
684:     PetscCall(MatSetValuesLocal(jac, 1, &i, 1 + user->itot[i], col, value, INSERT_VALUES));
685:   }

687:   /*
688:      Assemble matrix, using the 2-step process:
689:        MatAssemblyBegin(), MatAssemblyEnd().
690:      Between these two calls, the pointer to vector data has been restored to
691:      demonstrate the use of overlapping communicationn with computation.
692:   */
693:   PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY));
694:   PetscCall(VecRestoreArray(localX, &x));
695:   PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY));

697:   /*
698:      Tell the matrix we will never add a new nonzero location to the
699:      matrix. If we do, it will generate an error.
700:   */
701:   PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
702:   /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */
703:   PetscFunctionReturn(PETSC_SUCCESS);
704: }

706: /*TEST

708:    build:
709:       requires: !complex

711:    test:
712:       nsize: 2
713:       args: -snes_monitor_short
714:       localrunfiles: options.inf adj.in

716:    test:
717:       suffix: 2
718:       nsize: 2
719:       args: -snes_monitor_short -fd_jacobian_coloring
720:       localrunfiles: options.inf adj.in

722: TEST*/