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