1 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> 13 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"; 32 33 /* ------------------------------------------------------------------------ 34 35 PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian. 36 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. 40 41 -----------------------------------------------------------------------*/ 42 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 */ 46 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; 64 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); 71 72 int main(int argc, char **argv) { 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 sprintf(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 PetscInt i, Nvlocal; 498 PetscInt *gloInd; 499 PetscScalar *x; 500 #if defined(UNUSED_VARIABLES) 501 PetscReal temp1, temp, hx, hy, hxdhy, hydhx, sc; 502 PetscInt Neglobal, Nvglobal, j, row; 503 PetscReal alpha, lambda; 504 505 Nvglobal = user->Nvglobal; 506 Neglobal = user->Neglobal; 507 lambda = user->non_lin_param; 508 alpha = user->lin_param; 509 #endif 510 511 Nvlocal = user->Nvlocal; 512 gloInd = user->gloInd; 513 514 /* 515 Get a pointer to vector data. 516 - For default PETSc vectors, VecGetArray() returns a pointer to 517 the data array. Otherwise, the routine is implementation dependent. 518 - You MUST call VecRestoreArray() when you no longer need access to 519 the array. 520 */ 521 PetscCall(VecGetArray(X, &x)); 522 523 /* 524 Compute initial guess over the locally owned part of the grid 525 */ 526 for (i = 0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i]; 527 528 /* 529 Restore vector 530 */ 531 PetscCall(VecRestoreArray(X, &x)); 532 return 0; 533 } 534 /* -------------------- Evaluate Function F(x) --------------------- */ 535 /* 536 FormFunction - Evaluates nonlinear function, F(x). 537 538 Input Parameters: 539 . snes - the SNES context 540 . X - input vector 541 . ptr - optional user-defined context, as set by SNESSetFunction() 542 543 Output Parameter: 544 . F - function vector 545 */ 546 PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *ptr) { 547 AppCtx *user = (AppCtx *)ptr; 548 PetscInt i, j, Nvlocal; 549 PetscReal alpha, lambda; 550 PetscScalar *x, *f; 551 VecScatter scatter; 552 Vec localX = user->localX; 553 #if defined(UNUSED_VARIABLES) 554 PetscScalar ut, ub, ul, ur, u, *g, sc, uyy, uxx; 555 PetscReal hx, hy, hxdhy, hydhx; 556 PetscReal two = 2.0, one = 1.0; 557 PetscInt Nvglobal, Neglobal, row; 558 PetscInt *gloInd; 559 560 Nvglobal = user->Nvglobal; 561 Neglobal = user->Neglobal; 562 gloInd = user->gloInd; 563 #endif 564 565 Nvlocal = user->Nvlocal; 566 lambda = user->non_lin_param; 567 alpha = user->lin_param; 568 scatter = user->scatter; 569 570 /* 571 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 572 described in the beginning of this code 573 574 First scatter the distributed vector X into local vector localX (that includes 575 values for ghost nodes. If we wish,we can put some other work between 576 VecScatterBegin() and VecScatterEnd() to overlap the communication with 577 computation. 578 */ 579 PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 580 PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 581 582 /* 583 Get pointers to vector data 584 */ 585 PetscCall(VecGetArray(localX, &x)); 586 PetscCall(VecGetArray(F, &f)); 587 588 /* 589 Now compute the f(x). As mentioned earlier, the computed Laplacian is just an 590 approximate one chosen for illustrative purpose only. Another point to notice 591 is that this is a local (completly parallel) calculation. In practical application 592 codes, function calculation time is a dominat portion of the overall execution time. 593 */ 594 for (i = 0; i < Nvlocal; i++) { 595 f[i] = (user->itot[i] - alpha) * x[i] - lambda * x[i] * x[i]; 596 for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]]; 597 } 598 599 /* 600 Restore vectors 601 */ 602 PetscCall(VecRestoreArray(localX, &x)); 603 PetscCall(VecRestoreArray(F, &f)); 604 /*PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD));*/ 605 606 return 0; 607 } 608 609 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 610 /* 611 FormJacobian - Evaluates Jacobian matrix. 612 613 Input Parameters: 614 . snes - the SNES context 615 . X - input vector 616 . ptr - optional user-defined context, as set by SNESSetJacobian() 617 618 Output Parameters: 619 . A - Jacobian matrix 620 . B - optionally different preconditioning matrix 621 . flag - flag indicating matrix structure 622 623 */ 624 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr) { 625 AppCtx *user = (AppCtx *)ptr; 626 PetscInt i, j, Nvlocal, col[50]; 627 PetscScalar alpha, lambda, value[50]; 628 Vec localX = user->localX; 629 VecScatter scatter; 630 PetscScalar *x; 631 #if defined(UNUSED_VARIABLES) 632 PetscScalar two = 2.0, one = 1.0; 633 PetscInt row, Nvglobal, Neglobal; 634 PetscInt *gloInd; 635 636 Nvglobal = user->Nvglobal; 637 Neglobal = user->Neglobal; 638 gloInd = user->gloInd; 639 #endif 640 641 /*printf("Entering into FormJacobian \n");*/ 642 Nvlocal = user->Nvlocal; 643 lambda = user->non_lin_param; 644 alpha = user->lin_param; 645 scatter = user->scatter; 646 647 /* 648 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 649 described in the beginning of this code 650 651 First scatter the distributed vector X into local vector localX (that includes 652 values for ghost nodes. If we wish, we can put some other work between 653 VecScatterBegin() and VecScatterEnd() to overlap the communication with 654 computation. 655 */ 656 PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 657 PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 658 659 /* 660 Get pointer to vector data 661 */ 662 PetscCall(VecGetArray(localX, &x)); 663 664 for (i = 0; i < Nvlocal; i++) { 665 col[0] = i; 666 value[0] = user->itot[i] - 2.0 * lambda * x[i] - alpha; 667 for (j = 0; j < user->itot[i]; j++) { 668 col[j + 1] = user->AdjM[i][j]; 669 value[j + 1] = -1.0; 670 } 671 672 /* 673 Set the matrix values in the local ordering. Note that in order to use this 674 feature we must call the routine MatSetLocalToGlobalMapping() after the 675 matrix has been created. 676 */ 677 PetscCall(MatSetValuesLocal(jac, 1, &i, 1 + user->itot[i], col, value, INSERT_VALUES)); 678 } 679 680 /* 681 Assemble matrix, using the 2-step process: 682 MatAssemblyBegin(), MatAssemblyEnd(). 683 Between these two calls, the pointer to vector data has been restored to 684 demonstrate the use of overlapping communicationn with computation. 685 */ 686 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 687 PetscCall(VecRestoreArray(localX, &x)); 688 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 689 690 /* 691 Tell the matrix we will never add a new nonzero location to the 692 matrix. If we do, it will generate an error. 693 */ 694 PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 695 /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */ 696 return 0; 697 } 698 699 /*TEST 700 701 build: 702 requires: !complex 703 704 test: 705 nsize: 2 706 args: -snes_monitor_short 707 localrunfiles: options.inf adj.in 708 709 test: 710 suffix: 2 711 nsize: 2 712 args: -snes_monitor_short -fd_jacobian_coloring 713 localrunfiles: options.inf adj.in 714 715 TEST*/ 716