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 { 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; 105 106 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 107 Initialize program 108 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 109 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)); 114 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."); 117 118 /* 119 Initialize problem parameters 120 */ 121 user.Nvglobal = 16; /* Global # of vertices */ 122 user.Neglobal = 18; /* Global # of elements */ 123 124 PetscCall(PetscOptionsGetInt(NULL, NULL, "-vert", &user.Nvglobal, NULL)); 125 PetscCall(PetscOptionsGetInt(NULL, NULL, "-elem", &user.Neglobal, NULL)); 126 127 user.non_lin_param = 0.06; 128 129 PetscCall(PetscOptionsGetReal(NULL, NULL, "-nl_par", &user.non_lin_param, NULL)); 130 131 user.lin_param = -1.0; 132 133 PetscCall(PetscOptionsGetReal(NULL, NULL, "-lin_par", &user.lin_param, NULL)); 134 135 user.Nvlocal = 0; 136 user.Nelocal = 0; 137 138 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 139 Read the mesh and partitioning information 140 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 141 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. 150 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"); 159 160 /* 161 Each processor writes to the file output.<rank> where rank is the 162 processor's rank. 163 */ 164 sprintf(part_name, "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])); 175 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)); 180 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))); 187 188 sscanf(str, form, &dtmp); 189 user.AdjM[user.Nvlocal][i] = dtmp; 190 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)); 198 199 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 200 Create different orderings 201 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 202 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)); 212 213 for (i = 0; i < user.Nvlocal; i++) pordering[i] = rstart + i; 214 215 /* 216 Create the AO object 217 */ 218 PetscCall(AOCreateBasic(MPI_COMM_WORLD, user.Nvlocal, user.gloInd, pordering, &ao)); 219 PetscCall(PetscFree(pordering)); 220 221 /* 222 Keep the global indices for later use 223 */ 224 PetscCall(PetscMalloc1(user.Nvlocal, &user.locInd)); 225 PetscCall(PetscMalloc1(Nvneighborstotal, &tmp)); 226 227 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 228 Demonstrate the use of AO functionality 229 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 230 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])); 234 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])); 243 244 tmp[j + jstart] = user.AdjM[i][j]; 245 } 246 jstart += user.itot[i]; 247 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n")); 248 } 249 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)); 256 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")); 260 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]; 266 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 } 272 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 286 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; 295 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 } 303 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 } 316 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")); 321 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; 333 334 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "%" PetscInt_FMT " ", user.AdjM[i][j])); 335 } 336 PetscCall(PetscFPrintf(PETSC_COMM_SELF, fptr1, "\n")); 337 } 338 339 N = user.Nvglobal; 340 341 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 342 Create vector and matrix data structures 343 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 344 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)); 354 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)); 364 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)); 374 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)); 380 381 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 382 Create nonlinear solver context 383 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 384 385 PetscCall(SNESCreate(MPI_COMM_WORLD, &snes)); 386 PetscCall(SNESSetType(snes, type)); 387 388 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 389 Set routines for function and Jacobian evaluation 390 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 391 PetscCall(SNESSetFunction(snes, r, FormFunction, (void *)&user)); 392 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; 399 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 } 416 417 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 418 Customize nonlinear solver; set runtime options 419 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 420 421 PetscCall(SNESSetFromOptions(snes)); 422 423 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 424 Evaluate initial guess; then solve nonlinear system 425 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 426 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)); 434 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)); 441 442 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 443 Now solve the nonlinear system 444 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 445 446 PetscCall(SNESSolve(snes, NULL, x)); 447 PetscCall(SNESGetIterationNumber(snes, &its)); 448 PetscCall(SNESGetNonlinearStepFailures(snes, &nfails)); 449 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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 455 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)); 462 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 ----------------- */ 486 487 /* 488 FormInitialGuess - Forms initial approximation. 489 490 Input Parameters: 491 user - user-defined application context 492 X - vector 493 494 Output Parameter: 495 X - vector 496 */ 497 PetscErrorCode FormInitialGuess(AppCtx *user, Vec X) 498 { 499 PetscInt i, 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; 506 507 Nvglobal = user->Nvglobal; 508 Neglobal = user->Neglobal; 509 lambda = user->non_lin_param; 510 alpha = user->lin_param; 511 #endif 512 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 (i = 0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i]; 529 530 /* 531 Restore vector 532 */ 533 PetscCall(VecRestoreArray(X, &x)); 534 return 0; 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 i, j, 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 Nvlocal = user->Nvlocal; 569 lambda = user->non_lin_param; 570 alpha = user->lin_param; 571 scatter = user->scatter; 572 573 /* 574 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 575 described in the beginning of this code 576 577 First scatter the distributed vector X into local vector localX (that includes 578 values for ghost nodes. If we wish,we can put some other work between 579 VecScatterBegin() and VecScatterEnd() to overlap the communication with 580 computation. 581 */ 582 PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 583 PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 584 585 /* 586 Get pointers to vector data 587 */ 588 PetscCall(VecGetArray(localX, &x)); 589 PetscCall(VecGetArray(F, &f)); 590 591 /* 592 Now compute the f(x). As mentioned earlier, the computed Laplacian is just an 593 approximate one chosen for illustrative purpose only. Another point to notice 594 is that this is a local (completly parallel) calculation. In practical application 595 codes, function calculation time is a dominat portion of the overall execution time. 596 */ 597 for (i = 0; i < Nvlocal; i++) { 598 f[i] = (user->itot[i] - alpha) * x[i] - lambda * x[i] * x[i]; 599 for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]]; 600 } 601 602 /* 603 Restore vectors 604 */ 605 PetscCall(VecRestoreArray(localX, &x)); 606 PetscCall(VecRestoreArray(F, &f)); 607 /* PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD)); */ 608 609 return 0; 610 } 611 612 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 613 /* 614 FormJacobian - Evaluates Jacobian matrix. 615 616 Input Parameters: 617 . snes - the SNES context 618 . X - input vector 619 . ptr - optional user-defined context, as set by SNESSetJacobian() 620 621 Output Parameters: 622 . A - Jacobian matrix 623 . B - optionally different preconditioning matrix 624 . flag - flag indicating matrix structure 625 626 */ 627 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat jac, void *ptr) 628 { 629 AppCtx *user = (AppCtx *)ptr; 630 PetscInt i, j, Nvlocal, col[50]; 631 PetscScalar alpha, lambda, value[50]; 632 Vec localX = user->localX; 633 VecScatter scatter; 634 PetscScalar *x; 635 #if defined(UNUSED_VARIABLES) 636 PetscScalar two = 2.0, one = 1.0; 637 PetscInt row, Nvglobal, Neglobal; 638 PetscInt *gloInd; 639 640 Nvglobal = user->Nvglobal; 641 Neglobal = user->Neglobal; 642 gloInd = user->gloInd; 643 #endif 644 645 /*printf("Entering into FormJacobian \n");*/ 646 Nvlocal = user->Nvlocal; 647 lambda = user->non_lin_param; 648 alpha = user->lin_param; 649 scatter = user->scatter; 650 651 /* 652 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 653 described in the beginning of this code 654 655 First scatter the distributed vector X into local vector localX (that includes 656 values for ghost nodes. If we wish, we can put some other work between 657 VecScatterBegin() and VecScatterEnd() to overlap the communication with 658 computation. 659 */ 660 PetscCall(VecScatterBegin(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 661 PetscCall(VecScatterEnd(scatter, X, localX, INSERT_VALUES, SCATTER_FORWARD)); 662 663 /* 664 Get pointer to vector data 665 */ 666 PetscCall(VecGetArray(localX, &x)); 667 668 for (i = 0; i < Nvlocal; i++) { 669 col[0] = i; 670 value[0] = user->itot[i] - 2.0 * lambda * x[i] - alpha; 671 for (j = 0; j < user->itot[i]; j++) { 672 col[j + 1] = user->AdjM[i][j]; 673 value[j + 1] = -1.0; 674 } 675 676 /* 677 Set the matrix values in the local ordering. Note that in order to use this 678 feature we must call the routine MatSetLocalToGlobalMapping() after the 679 matrix has been created. 680 */ 681 PetscCall(MatSetValuesLocal(jac, 1, &i, 1 + user->itot[i], col, value, INSERT_VALUES)); 682 } 683 684 /* 685 Assemble matrix, using the 2-step process: 686 MatAssemblyBegin(), MatAssemblyEnd(). 687 Between these two calls, the pointer to vector data has been restored to 688 demonstrate the use of overlapping communicationn with computation. 689 */ 690 PetscCall(MatAssemblyBegin(jac, MAT_FINAL_ASSEMBLY)); 691 PetscCall(VecRestoreArray(localX, &x)); 692 PetscCall(MatAssemblyEnd(jac, MAT_FINAL_ASSEMBLY)); 693 694 /* 695 Tell the matrix we will never add a new nonzero location to the 696 matrix. If we do, it will generate an error. 697 */ 698 PetscCall(MatSetOption(jac, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE)); 699 /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */ 700 return 0; 701 } 702 703 /*TEST 704 705 build: 706 requires: !complex 707 708 test: 709 nsize: 2 710 args: -snes_monitor_short 711 localrunfiles: options.inf adj.in 712 713 test: 714 suffix: 2 715 nsize: 2 716 args: -snes_monitor_short -fd_jacobian_coloring 717 localrunfiles: options.inf adj.in 718 719 TEST*/ 720