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 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, (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 } 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 */ 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 */ 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 610 PetscFunctionReturn(PETSC_SUCCESS); 611 } 612 613 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 614 /* 615 FormJacobian - Evaluates Jacobian matrix. 616 617 Input Parameters: 618 . snes - the SNES context 619 . X - input vector 620 . ptr - optional user-defined context, as set by SNESSetJacobian() 621 622 Output Parameters: 623 . A - Jacobian matrix 624 . B - optionally different preconditioning matrix 625 . flag - flag indicating matrix structure 626 627 */ 628 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr) 629 { 630 AppCtx *user = (AppCtx *)ptr; 631 PetscInt Nvlocal, col[50]; 632 PetscScalar alpha, lambda, value[50]; 633 Vec localX = user->localX; 634 VecScatter scatter; 635 PetscScalar *x; 636 #if defined(UNUSED_VARIABLES) 637 PetscScalar two = 2.0, one = 1.0; 638 PetscInt row, Nvglobal, Neglobal; 639 PetscInt *gloInd; 640 641 Nvglobal = user->Nvglobal; 642 Neglobal = user->Neglobal; 643 gloInd = user->gloInd; 644 #endif 645 646 PetscFunctionBeginUser; 647 /*printf("Entering into FormJacobian \n");*/ 648 Nvlocal = user->Nvlocal; 649 lambda = user->non_lin_param; 650 alpha = user->lin_param; 651 scatter = user->scatter; 652 653 /* 654 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 655 described in the beginning of this code 656 657 First scatter the distributed vector X into local vector localX (that includes 658 values for ghost nodes. If we wish, we can put some other work between 659 VecScatterBegin() and VecScatterEnd() to overlap the communication with 660 computation. 661 */ 662 PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 663 PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 664 665 /* 666 Get pointer to vector data 667 */ 668 PetscCall(VecGetArray(localX, &x)); 669 670 for (PetscInt i = 0; i < Nvlocal; i++) { 671 col[0] = i; 672 value[0] = user->itot[i] - 2.0 * lambda * x[i] - alpha; 673 for (PetscInt j = 0; j < user->itot[i]; j++) { 674 col[j + 1] = user->AdjM[i][j]; 675 value[j + 1] = -1.0; 676 } 677 678 /* 679 Set the matrix values in the local ordering. Note that in order to use this 680 feature we must call the routine MatSetLocalToGlobalMapping() after the 681 matrix has been created. 682 */ 683 PetscCall(MatSetValuesLocal(jac, 1, &i, 1 + user->itot[i], col, value, INSERT_VALUES)); 684 } 685 686 /* 687 Assemble matrix, using the 2-step process: 688 MatAssemblyBegin(), MatAssemblyEnd(). 689 Between these two calls, the pointer to vector data has been restored to 690 demonstrate the use of overlapping communicationn with computation. 691 */ 692 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 693 PetscCall(VecRestoreArray(localX, &x)); 694 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 695 696 /* 697 Tell the matrix we will never add a new nonzero location to the 698 matrix. If we do, it will generate an error. 699 */ 700 PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 701 /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */ 702 PetscFunctionReturn(PETSC_SUCCESS); 703 } 704 705 /*TEST 706 707 build: 708 requires: !complex 709 710 test: 711 nsize: 2 712 args: -snes_monitor_short 713 localrunfiles: options.inf adj.in 714 715 test: 716 suffix: 2 717 nsize: 2 718 args: -snes_monitor_short -fd_jacobian_coloring 719 localrunfiles: options.inf adj.in 720 721 TEST*/ 722