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

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

    639:   Nvglobal = user->Nvglobal;
    640:   Neglobal = user->Neglobal;
    641:   gloInd   = user->gloInd;
    642: #endif

    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;

    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:   PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));
    661:   PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD));

    663:   /*
    664:      Get pointer to vector data
    665:   */
    666:   PetscCall(VecGetArray(localX, &x));

    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:     }

    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:   }

    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));

    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: }

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