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 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);user.v2p[inode] = dtmp; 171 if (user.v2p[inode] == rank) { 172 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %" PetscInt_FMT " belongs to processor %" PetscInt_FMT "\n",inode,user.v2p[inode])); 173 174 user.gloInd[user.Nvlocal] = inode; 175 sscanf(str,"%*d %d",&dtmp); 176 nbrs = dtmp; 177 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"Number of neighbors for the vertex %" PetscInt_FMT " is %" PetscInt_FMT "\n",inode,nbrs)); 178 179 user.itot[user.Nvlocal] = nbrs; 180 Nvneighborstotal += nbrs; 181 for (i = 0; i < user.itot[user.Nvlocal]; i++) { 182 form[0]='\0'; 183 for (j=0; j < i+2; j++) { 184 PetscCall(PetscStrlcat(form,"%*d ",sizeof(form))); 185 } 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++) { 259 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1," %" PetscInt_FMT " ",user.locInd[i])); 260 } 261 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n")); 262 263 jstart = 0; 264 for (i=0; i < user.Nvlocal; i++) { 265 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %" PetscInt_FMT " are : ",user.locInd[i])); 266 for (j=0; j < user.itot[i]; j++) { 267 user.AdjM[i][j] = tmp[j+jstart]; 268 269 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"%" PetscInt_FMT " ",user.AdjM[i][j])); 270 } 271 jstart += user.itot[i]; 272 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n")); 273 } 274 275 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 276 Extract the ghost vertex information for each processor 277 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 278 /* 279 Next, we need to generate a list of vertices required for this processor 280 and a local numbering scheme for all vertices required on this processor. 281 vertices - integer array of all vertices needed on this processor in PETSc 282 global numbering; this list consists of first the "locally owned" 283 vertices followed by the ghost vertices. 284 verticesmask - integer array that for each global vertex lists its local 285 vertex number (in vertices) + 1. If the global vertex is not 286 represented on this processor, then the corresponding 287 entry in verticesmask is zero 288 289 Note: vertices and verticesmask are both Nvglobal in length; this may 290 sound terribly non-scalable, but in fact is not so bad for a reasonable 291 number of processors. Importantly, it allows us to use NO SEARCHING 292 in setting up the data structures. 293 */ 294 PetscCall(PetscMalloc1(user.Nvglobal,&vertices)); 295 PetscCall(PetscCalloc1(user.Nvglobal,&verticesmask)); 296 nvertices = 0; 297 298 /* 299 First load "owned vertices" into list 300 */ 301 for (i=0; i < user.Nvlocal; i++) { 302 vertices[nvertices++] = user.locInd[i]; 303 verticesmask[user.locInd[i]] = nvertices; 304 } 305 306 /* 307 Now load ghost vertices into list 308 */ 309 for (i=0; i < user.Nvlocal; i++) { 310 for (j=0; j < user.itot[i]; j++) { 311 nb = user.AdjM[i][j]; 312 if (!verticesmask[nb]) { 313 vertices[nvertices++] = nb; 314 verticesmask[nb] = nvertices; 315 } 316 } 317 } 318 319 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n")); 320 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n")); 321 for (i=0; i < nvertices; i++) { 322 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"%" PetscInt_FMT " ",vertices[i])); 323 } 324 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n")); 325 326 /* 327 Map the vertices listed in the neighbors to the local numbering from 328 the global ordering that they contained initially. 329 */ 330 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n")); 331 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"After mapping neighbors in the local contiguous ordering\n")); 332 for (i=0; i<user.Nvlocal; i++) { 333 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %" PetscInt_FMT " are :\n",i)); 334 for (j = 0; j < user.itot[i]; j++) { 335 nb = user.AdjM[i][j]; 336 user.AdjM[i][j] = verticesmask[nb] - 1; 337 338 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"%" PetscInt_FMT " ",user.AdjM[i][j])); 339 } 340 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n")); 341 } 342 343 N = user.Nvglobal; 344 345 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 346 Create vector and matrix data structures 347 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 348 349 /* 350 Create vector data structures 351 */ 352 PetscCall(VecCreate(MPI_COMM_WORLD,&x)); 353 PetscCall(VecSetSizes(x,user.Nvlocal,N)); 354 PetscCall(VecSetFromOptions(x)); 355 PetscCall(VecDuplicate(x,&r)); 356 PetscCall(VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX)); 357 PetscCall(VecDuplicate(user.localX,&user.localF)); 358 359 /* 360 Create the scatter between the global representation and the 361 local representation 362 */ 363 PetscCall(ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal)); 364 PetscCall(ISCreateBlock(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isglobal)); 365 PetscCall(VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter)); 366 PetscCall(ISDestroy(&isglobal)); 367 PetscCall(ISDestroy(&islocal)); 368 369 /* 370 Create matrix data structure; Just to keep the example simple, we have not done any 371 preallocation of memory for the matrix. In real application code with big matrices, 372 preallocation should always be done to expedite the matrix creation. 373 */ 374 PetscCall(MatCreate(MPI_COMM_WORLD,&Jac)); 375 PetscCall(MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N)); 376 PetscCall(MatSetFromOptions(Jac)); 377 PetscCall(MatSetUp(Jac)); 378 379 /* 380 The following routine allows us to set the matrix values in local ordering 381 */ 382 PetscCall(ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isl2g)); 383 PetscCall(MatSetLocalToGlobalMapping(Jac,isl2g,isl2g)); 384 385 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 386 Create nonlinear solver context 387 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 388 389 PetscCall(SNESCreate(MPI_COMM_WORLD,&snes)); 390 PetscCall(SNESSetType(snes,type)); 391 392 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 393 Set routines for function and Jacobian evaluation 394 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 395 PetscCall(SNESSetFunction(snes,r,FormFunction,(void*)&user)); 396 397 PetscCall(PetscOptionsGetBool(NULL,NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,0)); 398 if (!fd_jacobian_coloring) { 399 PetscCall(SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user)); 400 } else { /* Use matfdcoloring */ 401 ISColoring iscoloring; 402 MatColoring mc; 403 404 /* Get the data structure of Jac */ 405 PetscCall(FormJacobian(snes,x,Jac,Jac,&user)); 406 /* Create coloring context */ 407 PetscCall(MatColoringCreate(Jac,&mc)); 408 PetscCall(MatColoringSetType(mc,MATCOLORINGSL)); 409 PetscCall(MatColoringSetFromOptions(mc)); 410 PetscCall(MatColoringApply(mc,&iscoloring)); 411 PetscCall(MatColoringDestroy(&mc)); 412 PetscCall(MatFDColoringCreate(Jac,iscoloring,&matfdcoloring)); 413 PetscCall(MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user)); 414 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 415 PetscCall(MatFDColoringSetUp(Jac,iscoloring,matfdcoloring)); 416 /* PetscCall(MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD)); */ 417 PetscCall(SNESSetJacobian(snes,Jac,Jac,SNESComputeJacobianDefaultColor,matfdcoloring)); 418 PetscCall(ISColoringDestroy(&iscoloring)); 419 } 420 421 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 422 Customize nonlinear solver; set runtime options 423 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 424 425 PetscCall(SNESSetFromOptions(snes)); 426 427 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 428 Evaluate initial guess; then solve nonlinear system 429 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 430 431 /* 432 Note: The user should initialize the vector, x, with the initial guess 433 for the nonlinear solver prior to calling SNESSolve(). In particular, 434 to employ an initial guess of zero, the user should explicitly set 435 this vector to zero by calling VecSet(). 436 */ 437 PetscCall(FormInitialGuess(&user,x)); 438 439 /* 440 Print the initial guess 441 */ 442 PetscCall(VecGetArray(x,&xx)); 443 for (inode = 0; inode < user.Nvlocal; inode++) { 444 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"Initial Solution at node %" PetscInt_FMT " is %f \n",inode,(double)PetscRealPart(xx[inode]))); 445 } 446 PetscCall(VecRestoreArray(x,&xx)); 447 448 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 449 Now solve the nonlinear system 450 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 451 452 PetscCall(SNESSolve(snes,NULL,x)); 453 PetscCall(SNESGetIterationNumber(snes,&its)); 454 PetscCall(SNESGetNonlinearStepFailures(snes,&nfails)); 455 456 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 457 Print the output : solution vector and other information 458 Each processor writes to the file output.<rank> where rank is the 459 processor's rank. 460 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 461 462 PetscCall(VecGetArray(x,&xx)); 463 for (inode = 0; inode < user.Nvlocal; inode++) { 464 PetscCall(PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %" PetscInt_FMT " is %f \n",inode,(double)PetscRealPart(xx[inode]))); 465 } 466 PetscCall(VecRestoreArray(x,&xx)); 467 fclose(fptr1); 468 PetscCall(PetscPrintf(MPI_COMM_WORLD,"number of SNES iterations = %" PetscInt_FMT ", ",its)); 469 PetscCall(PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %" PetscInt_FMT "\n",nfails)); 470 471 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 472 Free work space. All PETSc objects should be destroyed when they 473 are no longer needed. 474 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 475 PetscCall(PetscFree(user.gloInd)); 476 PetscCall(PetscFree(user.locInd)); 477 PetscCall(PetscFree(vertices)); 478 PetscCall(PetscFree(verticesmask)); 479 PetscCall(PetscFree(tmp)); 480 PetscCall(VecScatterDestroy(&user.scatter)); 481 PetscCall(ISLocalToGlobalMappingDestroy(&isl2g)); 482 PetscCall(VecDestroy(&x)); 483 PetscCall(VecDestroy(&r)); 484 PetscCall(VecDestroy(&user.localX)); 485 PetscCall(VecDestroy(&user.localF)); 486 PetscCall(SNESDestroy(&snes)); 487 PetscCall(MatDestroy(&Jac)); 488 /* PetscCall(PetscDrawDestroy(draw));*/ 489 if (fd_jacobian_coloring) PetscCall(MatFDColoringDestroy(&matfdcoloring)); 490 PetscCall(PetscFinalize()); 491 return 0; 492 } 493 /* -------------------- Form initial approximation ----------------- */ 494 495 /* 496 FormInitialGuess - Forms initial approximation. 497 498 Input Parameters: 499 user - user-defined application context 500 X - vector 501 502 Output Parameter: 503 X - vector 504 */ 505 PetscErrorCode FormInitialGuess(AppCtx *user,Vec X) 506 { 507 PetscInt i,Nvlocal; 508 PetscInt *gloInd; 509 PetscScalar *x; 510 #if defined(UNUSED_VARIABLES) 511 PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc; 512 PetscInt Neglobal,Nvglobal,j,row; 513 PetscReal alpha,lambda; 514 515 Nvglobal = user->Nvglobal; 516 Neglobal = user->Neglobal; 517 lambda = user->non_lin_param; 518 alpha = user->lin_param; 519 #endif 520 521 Nvlocal = user->Nvlocal; 522 gloInd = user->gloInd; 523 524 /* 525 Get a pointer to vector data. 526 - For default PETSc vectors, VecGetArray() returns a pointer to 527 the data array. Otherwise, the routine is implementation dependent. 528 - You MUST call VecRestoreArray() when you no longer need access to 529 the array. 530 */ 531 PetscCall(VecGetArray(X,&x)); 532 533 /* 534 Compute initial guess over the locally owned part of the grid 535 */ 536 for (i=0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i]; 537 538 /* 539 Restore vector 540 */ 541 PetscCall(VecRestoreArray(X,&x)); 542 return 0; 543 } 544 /* -------------------- Evaluate Function F(x) --------------------- */ 545 /* 546 FormFunction - Evaluates nonlinear function, F(x). 547 548 Input Parameters: 549 . snes - the SNES context 550 . X - input vector 551 . ptr - optional user-defined context, as set by SNESSetFunction() 552 553 Output Parameter: 554 . F - function vector 555 */ 556 PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr) 557 { 558 AppCtx *user = (AppCtx*)ptr; 559 PetscInt i,j,Nvlocal; 560 PetscReal alpha,lambda; 561 PetscScalar *x,*f; 562 VecScatter scatter; 563 Vec localX = user->localX; 564 #if defined(UNUSED_VARIABLES) 565 PetscScalar ut,ub,ul,ur,u,*g,sc,uyy,uxx; 566 PetscReal hx,hy,hxdhy,hydhx; 567 PetscReal two = 2.0,one = 1.0; 568 PetscInt Nvglobal,Neglobal,row; 569 PetscInt *gloInd; 570 571 Nvglobal = user->Nvglobal; 572 Neglobal = user->Neglobal; 573 gloInd = user->gloInd; 574 #endif 575 576 Nvlocal = user->Nvlocal; 577 lambda = user->non_lin_param; 578 alpha = user->lin_param; 579 scatter = user->scatter; 580 581 /* 582 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 583 described in the beginning of this code 584 585 First scatter the distributed vector X into local vector localX (that includes 586 values for ghost nodes. If we wish,we can put some other work between 587 VecScatterBegin() and VecScatterEnd() to overlap the communication with 588 computation. 589 */ 590 PetscCall(VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD)); 591 PetscCall(VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD)); 592 593 /* 594 Get pointers to vector data 595 */ 596 PetscCall(VecGetArray(localX,&x)); 597 PetscCall(VecGetArray(F,&f)); 598 599 /* 600 Now compute the f(x). As mentioned earlier, the computed Laplacian is just an 601 approximate one chosen for illustrative purpose only. Another point to notice 602 is that this is a local (completly parallel) calculation. In practical application 603 codes, function calculation time is a dominat portion of the overall execution time. 604 */ 605 for (i=0; i < Nvlocal; i++) { 606 f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i]; 607 for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]]; 608 } 609 610 /* 611 Restore vectors 612 */ 613 PetscCall(VecRestoreArray(localX,&x)); 614 PetscCall(VecRestoreArray(F,&f)); 615 /*PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD));*/ 616 617 return 0; 618 } 619 620 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 621 /* 622 FormJacobian - Evaluates Jacobian matrix. 623 624 Input Parameters: 625 . snes - the SNES context 626 . X - input vector 627 . ptr - optional user-defined context, as set by SNESSetJacobian() 628 629 Output Parameters: 630 . A - Jacobian matrix 631 . B - optionally different preconditioning matrix 632 . flag - flag indicating matrix structure 633 634 */ 635 PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr) 636 { 637 AppCtx *user = (AppCtx*)ptr; 638 PetscInt i,j,Nvlocal,col[50]; 639 PetscScalar alpha,lambda,value[50]; 640 Vec localX = user->localX; 641 VecScatter scatter; 642 PetscScalar *x; 643 #if defined(UNUSED_VARIABLES) 644 PetscScalar two = 2.0,one = 1.0; 645 PetscInt row,Nvglobal,Neglobal; 646 PetscInt *gloInd; 647 648 Nvglobal = user->Nvglobal; 649 Neglobal = user->Neglobal; 650 gloInd = user->gloInd; 651 #endif 652 653 /*printf("Entering into FormJacobian \n");*/ 654 Nvlocal = user->Nvlocal; 655 lambda = user->non_lin_param; 656 alpha = user->lin_param; 657 scatter = user->scatter; 658 659 /* 660 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 661 described in the beginning of this code 662 663 First scatter the distributed vector X into local vector localX (that includes 664 values for ghost nodes. If we wish, we can put some other work between 665 VecScatterBegin() and VecScatterEnd() to overlap the communication with 666 computation. 667 */ 668 PetscCall(VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD)); 669 PetscCall(VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD)); 670 671 /* 672 Get pointer to vector data 673 */ 674 PetscCall(VecGetArray(localX,&x)); 675 676 for (i=0; i < Nvlocal; i++) { 677 col[0] = i; 678 value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha; 679 for (j = 0; j < user->itot[i]; j++) { 680 col[j+1] = user->AdjM[i][j]; 681 value[j+1] = -1.0; 682 } 683 684 /* 685 Set the matrix values in the local ordering. Note that in order to use this 686 feature we must call the routine MatSetLocalToGlobalMapping() after the 687 matrix has been created. 688 */ 689 PetscCall(MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES)); 690 } 691 692 /* 693 Assemble matrix, using the 2-step process: 694 MatAssemblyBegin(), MatAssemblyEnd(). 695 Between these two calls, the pointer to vector data has been restored to 696 demonstrate the use of overlapping communicationn with computation. 697 */ 698 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 699 PetscCall(VecRestoreArray(localX,&x)); 700 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 701 702 /* 703 Tell the matrix we will never add a new nonzero location to the 704 matrix. If we do, it will generate an error. 705 */ 706 PetscCall(MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE)); 707 /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */ 708 return 0; 709 } 710 711 /*TEST 712 713 build: 714 requires: !complex 715 716 test: 717 nsize: 2 718 args: -snes_monitor_short 719 localrunfiles: options.inf adj.in 720 721 test: 722 suffix: 2 723 nsize: 2 724 args: -snes_monitor_short -fd_jacobian_coloring 725 localrunfiles: options.inf adj.in 726 727 TEST*/ 728