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 /*T 34 Concepts: SNES^unstructured grid 35 Concepts: AO^application to PETSc ordering or vice versa; 36 Concepts: VecScatter^using vector scatter operations; 37 Processors: n 38 T*/ 39 40 /* ------------------------------------------------------------------------ 41 42 PDE Solved : L(u) + lambda*u*u + alpha*u = 0 where L(u) is the Laplacian. 43 44 The Laplacian is approximated in the following way: each edge is given a weight 45 of one meaning that the diagonal term will have the weight equal to the degree 46 of a node. The off diagonal terms will get a weight of -1. 47 48 -----------------------------------------------------------------------*/ 49 50 #define MAX_ELEM 500 /* Maximum number of elements */ 51 #define MAX_VERT 100 /* Maximum number of vertices */ 52 #define MAX_VERT_ELEM 3 /* Vertices per element */ 53 54 /* 55 Application-defined context for problem specific data 56 */ 57 typedef struct { 58 PetscInt Nvglobal,Nvlocal; /* global and local number of vertices */ 59 PetscInt Neglobal,Nelocal; /* global and local number of vertices */ 60 PetscInt AdjM[MAX_VERT][50]; /* adjacency list of a vertex */ 61 PetscInt itot[MAX_VERT]; /* total number of neighbors for a vertex */ 62 PetscInt icv[MAX_ELEM][MAX_VERT_ELEM]; /* vertices belonging to an element */ 63 PetscInt v2p[MAX_VERT]; /* processor number for a vertex */ 64 PetscInt *locInd,*gloInd; /* local and global orderings for a node */ 65 Vec localX,localF; /* local solution (u) and f(u) vectors */ 66 PetscReal non_lin_param; /* nonlinear parameter for the PDE */ 67 PetscReal lin_param; /* linear parameter for the PDE */ 68 VecScatter scatter; /* scatter context for the local and 69 distributed vectors */ 70 } AppCtx; 71 72 /* 73 User-defined routines 74 */ 75 PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 76 PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 77 PetscErrorCode FormInitialGuess(AppCtx*,Vec); 78 79 int main(int argc,char **argv) 80 { 81 SNES snes; /* SNES context */ 82 SNESType type = SNESNEWTONLS; /* default nonlinear solution method */ 83 Vec x,r; /* solution, residual vectors */ 84 Mat Jac; /* Jacobian matrix */ 85 AppCtx user; /* user-defined application context */ 86 AO ao; /* Application Ordering object */ 87 IS isglobal,islocal; /* global and local index sets */ 88 PetscMPIInt rank,size; /* rank of a process, number of processors */ 89 PetscInt rstart; /* starting index of PETSc ordering for a processor */ 90 PetscInt nfails; /* number of unsuccessful Newton steps */ 91 PetscInt bs = 1; /* block size for multicomponent systems */ 92 PetscInt nvertices; /* number of local plus ghost nodes of a processor */ 93 PetscInt *pordering; /* PETSc ordering */ 94 PetscInt *vertices; /* list of all vertices (incl. ghost ones) on a processor */ 95 PetscInt *verticesmask; 96 PetscInt *tmp; 97 PetscInt i,j,jstart,inode,nb,nbrs,Nvneighborstotal = 0; 98 PetscErrorCode ierr; 99 PetscInt its,N; 100 PetscScalar *xx; 101 char str[256],form[256],part_name[256]; 102 FILE *fptr,*fptr1; 103 ISLocalToGlobalMapping isl2g; 104 int dtmp; 105 #if defined(UNUSED_VARIABLES) 106 PetscDraw draw; /* drawing context */ 107 PetscScalar *ff,*gg; 108 PetscReal tiny = 1.0e-10,zero = 0.0,one = 1.0,big = 1.0e+10; 109 PetscInt *tmp1,*tmp2; 110 #endif 111 MatFDColoring matfdcoloring = 0; 112 PetscBool fd_jacobian_coloring = PETSC_FALSE; 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Initialize program 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 118 ierr = PetscInitialize(&argc,&argv,"options.inf",help);if (ierr) return ierr; 119 ierr = MPI_Comm_rank(MPI_COMM_WORLD,&rank);CHKERRMPI(ierr); 120 ierr = MPI_Comm_size(MPI_COMM_WORLD,&size);CHKERRMPI(ierr); 121 122 /* The current input file options.inf is for 2 proc run only */ 123 if (size != 2) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"This example currently runs on 2 procs only."); 124 125 /* 126 Initialize problem parameters 127 */ 128 user.Nvglobal = 16; /*Global # of vertices */ 129 user.Neglobal = 18; /*Global # of elements */ 130 131 ierr = PetscOptionsGetInt(NULL,NULL,"-vert",&user.Nvglobal,NULL);CHKERRQ(ierr); 132 ierr = PetscOptionsGetInt(NULL,NULL,"-elem",&user.Neglobal,NULL);CHKERRQ(ierr); 133 134 user.non_lin_param = 0.06; 135 136 ierr = PetscOptionsGetReal(NULL,NULL,"-nl_par",&user.non_lin_param,NULL);CHKERRQ(ierr); 137 138 user.lin_param = -1.0; 139 140 ierr = PetscOptionsGetReal(NULL,NULL,"-lin_par",&user.lin_param,NULL);CHKERRQ(ierr); 141 142 user.Nvlocal = 0; 143 user.Nelocal = 0; 144 145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 Read the mesh and partitioning information 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 149 /* 150 Read the mesh and partitioning information from 'adj.in'. 151 The file format is as follows. 152 For each line the first entry is the processor rank where the 153 current node belongs. The second entry is the number of 154 neighbors of a node. The rest of the line is the adjacency 155 list of a node. Currently this file is set up to work on two 156 processors. 157 158 This is not a very good example of reading input. In the future, 159 we will put an example that shows the style that should be 160 used in a real application, where partitioning will be done 161 dynamically by calling partitioning routines (at present, we have 162 a ready interface to ParMeTiS). 163 */ 164 fptr = fopen("adj.in","r"); 165 if (!fptr) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Could not open adj.in"); 166 167 /* 168 Each processor writes to the file output.<rank> where rank is the 169 processor's rank. 170 */ 171 sprintf(part_name,"output.%d",rank); 172 fptr1 = fopen(part_name,"w"); 173 if (!fptr1) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Could no open output file"); 174 ierr = PetscMalloc1(user.Nvglobal,&user.gloInd);CHKERRQ(ierr); 175 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Rank is %d\n",rank);CHKERRQ(ierr); 176 for (inode = 0; inode < user.Nvglobal; inode++) { 177 if (!fgets(str,256,fptr)) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"fgets read failed"); 178 sscanf(str,"%d",&dtmp);user.v2p[inode] = dtmp; 179 if (user.v2p[inode] == rank) { 180 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Node %D belongs to processor %D\n",inode,user.v2p[inode]);CHKERRQ(ierr); 181 182 user.gloInd[user.Nvlocal] = inode; 183 sscanf(str,"%*d %d",&dtmp); 184 nbrs = dtmp; 185 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Number of neighbors for the vertex %D is %D\n",inode,nbrs);CHKERRQ(ierr); 186 187 user.itot[user.Nvlocal] = nbrs; 188 Nvneighborstotal += nbrs; 189 for (i = 0; i < user.itot[user.Nvlocal]; i++) { 190 form[0]='\0'; 191 for (j=0; j < i+2; j++) { 192 ierr = PetscStrlcat(form,"%*d ",sizeof(form));CHKERRQ(ierr); 193 } 194 ierr = PetscStrlcat(form,"%d",sizeof(form));CHKERRQ(ierr); 195 196 sscanf(str,form,&dtmp); 197 user.AdjM[user.Nvlocal][i] = dtmp; 198 199 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[user.Nvlocal][i]);CHKERRQ(ierr); 200 } 201 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 202 user.Nvlocal++; 203 } 204 } 205 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Total # of Local Vertices is %D \n",user.Nvlocal);CHKERRQ(ierr); 206 207 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 208 Create different orderings 209 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 210 211 /* 212 Create the local ordering list for vertices. First a list using the PETSc global 213 ordering is created. Then we use the AO object to get the PETSc-to-application and 214 application-to-PETSc mappings. Each vertex also gets a local index (stored in the 215 locInd array). 216 */ 217 ierr = MPI_Scan(&user.Nvlocal,&rstart,1,MPIU_INT,MPI_SUM,PETSC_COMM_WORLD);CHKERRMPI(ierr); 218 rstart -= user.Nvlocal; 219 ierr = PetscMalloc1(user.Nvlocal,&pordering);CHKERRQ(ierr); 220 221 for (i=0; i < user.Nvlocal; i++) pordering[i] = rstart + i; 222 223 /* 224 Create the AO object 225 */ 226 ierr = AOCreateBasic(MPI_COMM_WORLD,user.Nvlocal,user.gloInd,pordering,&ao);CHKERRQ(ierr); 227 ierr = PetscFree(pordering);CHKERRQ(ierr); 228 229 /* 230 Keep the global indices for later use 231 */ 232 ierr = PetscMalloc1(user.Nvlocal,&user.locInd);CHKERRQ(ierr); 233 ierr = PetscMalloc1(Nvneighborstotal,&tmp);CHKERRQ(ierr); 234 235 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 236 Demonstrate the use of AO functionality 237 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 238 239 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Before AOApplicationToPetsc, local indices are : \n");CHKERRQ(ierr); 240 for (i=0; i < user.Nvlocal; i++) { 241 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.gloInd[i]);CHKERRQ(ierr); 242 243 user.locInd[i] = user.gloInd[i]; 244 } 245 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 246 jstart = 0; 247 for (i=0; i < user.Nvlocal; i++) { 248 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.gloInd[i]);CHKERRQ(ierr); 249 for (j=0; j < user.itot[i]; j++) { 250 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);CHKERRQ(ierr); 251 252 tmp[j + jstart] = user.AdjM[i][j]; 253 } 254 jstart += user.itot[i]; 255 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 256 } 257 258 /* 259 Now map the vlocal and neighbor lists to the PETSc ordering 260 */ 261 ierr = AOApplicationToPetsc(ao,user.Nvlocal,user.locInd);CHKERRQ(ierr); 262 ierr = AOApplicationToPetsc(ao,Nvneighborstotal,tmp);CHKERRQ(ierr); 263 ierr = AODestroy(&ao);CHKERRQ(ierr); 264 265 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"After AOApplicationToPetsc, local indices are : \n");CHKERRQ(ierr); 266 for (i=0; i < user.Nvlocal; i++) { 267 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1," %D ",user.locInd[i]);CHKERRQ(ierr); 268 } 269 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 270 271 jstart = 0; 272 for (i=0; i < user.Nvlocal; i++) { 273 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are : ",user.locInd[i]);CHKERRQ(ierr); 274 for (j=0; j < user.itot[i]; j++) { 275 user.AdjM[i][j] = tmp[j+jstart]; 276 277 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);CHKERRQ(ierr); 278 } 279 jstart += user.itot[i]; 280 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 281 } 282 283 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 284 Extract the ghost vertex information for each processor 285 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 286 /* 287 Next, we need to generate a list of vertices required for this processor 288 and a local numbering scheme for all vertices required on this processor. 289 vertices - integer array of all vertices needed on this processor in PETSc 290 global numbering; this list consists of first the "locally owned" 291 vertices followed by the ghost vertices. 292 verticesmask - integer array that for each global vertex lists its local 293 vertex number (in vertices) + 1. If the global vertex is not 294 represented on this processor, then the corresponding 295 entry in verticesmask is zero 296 297 Note: vertices and verticesmask are both Nvglobal in length; this may 298 sound terribly non-scalable, but in fact is not so bad for a reasonable 299 number of processors. Importantly, it allows us to use NO SEARCHING 300 in setting up the data structures. 301 */ 302 ierr = PetscMalloc1(user.Nvglobal,&vertices);CHKERRQ(ierr); 303 ierr = PetscCalloc1(user.Nvglobal,&verticesmask);CHKERRQ(ierr); 304 nvertices = 0; 305 306 /* 307 First load "owned vertices" into list 308 */ 309 for (i=0; i < user.Nvlocal; i++) { 310 vertices[nvertices++] = user.locInd[i]; 311 verticesmask[user.locInd[i]] = nvertices; 312 } 313 314 /* 315 Now load ghost vertices into list 316 */ 317 for (i=0; i < user.Nvlocal; i++) { 318 for (j=0; j < user.itot[i]; j++) { 319 nb = user.AdjM[i][j]; 320 if (!verticesmask[nb]) { 321 vertices[nvertices++] = nb; 322 verticesmask[nb] = nvertices; 323 } 324 } 325 } 326 327 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 328 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"The array vertices is :\n");CHKERRQ(ierr); 329 for (i=0; i < nvertices; i++) { 330 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",vertices[i]);CHKERRQ(ierr); 331 } 332 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 333 334 /* 335 Map the vertices listed in the neighbors to the local numbering from 336 the global ordering that they contained initially. 337 */ 338 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 339 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"After mapping neighbors in the local contiguous ordering\n");CHKERRQ(ierr); 340 for (i=0; i<user.Nvlocal; i++) { 341 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Neghbors of local vertex %D are :\n",i);CHKERRQ(ierr); 342 for (j = 0; j < user.itot[i]; j++) { 343 nb = user.AdjM[i][j]; 344 user.AdjM[i][j] = verticesmask[nb] - 1; 345 346 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"%D ",user.AdjM[i][j]);CHKERRQ(ierr); 347 } 348 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"\n");CHKERRQ(ierr); 349 } 350 351 N = user.Nvglobal; 352 353 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 354 Create vector and matrix data structures 355 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 356 357 /* 358 Create vector data structures 359 */ 360 ierr = VecCreate(MPI_COMM_WORLD,&x);CHKERRQ(ierr); 361 ierr = VecSetSizes(x,user.Nvlocal,N);CHKERRQ(ierr); 362 ierr = VecSetFromOptions(x);CHKERRQ(ierr); 363 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 364 ierr = VecCreateSeq(MPI_COMM_SELF,bs*nvertices,&user.localX);CHKERRQ(ierr); 365 ierr = VecDuplicate(user.localX,&user.localF);CHKERRQ(ierr); 366 367 /* 368 Create the scatter between the global representation and the 369 local representation 370 */ 371 ierr = ISCreateStride(MPI_COMM_SELF,bs*nvertices,0,1,&islocal);CHKERRQ(ierr); 372 ierr = ISCreateBlock(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isglobal);CHKERRQ(ierr); 373 ierr = VecScatterCreate(x,isglobal,user.localX,islocal,&user.scatter);CHKERRQ(ierr); 374 ierr = ISDestroy(&isglobal);CHKERRQ(ierr); 375 ierr = ISDestroy(&islocal);CHKERRQ(ierr); 376 377 /* 378 Create matrix data structure; Just to keep the example simple, we have not done any 379 preallocation of memory for the matrix. In real application code with big matrices, 380 preallocation should always be done to expedite the matrix creation. 381 */ 382 ierr = MatCreate(MPI_COMM_WORLD,&Jac);CHKERRQ(ierr); 383 ierr = MatSetSizes(Jac,PETSC_DECIDE,PETSC_DECIDE,N,N);CHKERRQ(ierr); 384 ierr = MatSetFromOptions(Jac);CHKERRQ(ierr); 385 ierr = MatSetUp(Jac);CHKERRQ(ierr); 386 387 /* 388 The following routine allows us to set the matrix values in local ordering 389 */ 390 ierr = ISLocalToGlobalMappingCreate(MPI_COMM_SELF,bs,nvertices,vertices,PETSC_COPY_VALUES,&isl2g);CHKERRQ(ierr); 391 ierr = MatSetLocalToGlobalMapping(Jac,isl2g,isl2g);CHKERRQ(ierr); 392 393 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 394 Create nonlinear solver context 395 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 396 397 ierr = SNESCreate(MPI_COMM_WORLD,&snes);CHKERRQ(ierr); 398 ierr = SNESSetType(snes,type);CHKERRQ(ierr); 399 400 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 401 Set routines for function and Jacobian evaluation 402 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 403 ierr = SNESSetFunction(snes,r,FormFunction,(void*)&user);CHKERRQ(ierr); 404 405 ierr = PetscOptionsGetBool(NULL,NULL,"-fd_jacobian_coloring",&fd_jacobian_coloring,0);CHKERRQ(ierr); 406 if (!fd_jacobian_coloring) { 407 ierr = SNESSetJacobian(snes,Jac,Jac,FormJacobian,(void*)&user);CHKERRQ(ierr); 408 } else { /* Use matfdcoloring */ 409 ISColoring iscoloring; 410 MatColoring mc; 411 412 /* Get the data structure of Jac */ 413 ierr = FormJacobian(snes,x,Jac,Jac,&user);CHKERRQ(ierr); 414 /* Create coloring context */ 415 ierr = MatColoringCreate(Jac,&mc);CHKERRQ(ierr); 416 ierr = MatColoringSetType(mc,MATCOLORINGSL);CHKERRQ(ierr); 417 ierr = MatColoringSetFromOptions(mc);CHKERRQ(ierr); 418 ierr = MatColoringApply(mc,&iscoloring);CHKERRQ(ierr); 419 ierr = MatColoringDestroy(&mc);CHKERRQ(ierr); 420 ierr = MatFDColoringCreate(Jac,iscoloring,&matfdcoloring);CHKERRQ(ierr); 421 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);CHKERRQ(ierr); 422 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 423 ierr = MatFDColoringSetUp(Jac,iscoloring,matfdcoloring);CHKERRQ(ierr); 424 /* ierr = MatFDColoringView(matfdcoloring,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */ 425 ierr = SNESSetJacobian(snes,Jac,Jac,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr); 426 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 427 } 428 429 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 430 Customize nonlinear solver; set runtime options 431 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 432 433 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 434 435 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 436 Evaluate initial guess; then solve nonlinear system 437 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 438 439 /* 440 Note: The user should initialize the vector, x, with the initial guess 441 for the nonlinear solver prior to calling SNESSolve(). In particular, 442 to employ an initial guess of zero, the user should explicitly set 443 this vector to zero by calling VecSet(). 444 */ 445 ierr = FormInitialGuess(&user,x);CHKERRQ(ierr); 446 447 /* 448 Print the initial guess 449 */ 450 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 451 for (inode = 0; inode < user.Nvlocal; inode++) { 452 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Initial Solution at node %D is %f \n",inode,xx[inode]);CHKERRQ(ierr); 453 } 454 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 455 456 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 457 Now solve the nonlinear system 458 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 459 460 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); 461 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 462 ierr = SNESGetNonlinearStepFailures(snes,&nfails);CHKERRQ(ierr); 463 464 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 465 Print the output : solution vector and other information 466 Each processor writes to the file output.<rank> where rank is the 467 processor's rank. 468 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 469 470 ierr = VecGetArray(x,&xx);CHKERRQ(ierr); 471 for (inode = 0; inode < user.Nvlocal; inode++) { 472 ierr = PetscFPrintf(PETSC_COMM_SELF,fptr1,"Solution at node %D is %f \n",inode,xx[inode]);CHKERRQ(ierr); 473 } 474 ierr = VecRestoreArray(x,&xx);CHKERRQ(ierr); 475 fclose(fptr1); 476 ierr = PetscPrintf(MPI_COMM_WORLD,"number of SNES iterations = %D, ",its);CHKERRQ(ierr); 477 ierr = PetscPrintf(MPI_COMM_WORLD,"number of unsuccessful steps = %D\n",nfails);CHKERRQ(ierr); 478 479 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 480 Free work space. All PETSc objects should be destroyed when they 481 are no longer needed. 482 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 483 ierr = PetscFree(user.gloInd);CHKERRQ(ierr); 484 ierr = PetscFree(user.locInd);CHKERRQ(ierr); 485 ierr = PetscFree(vertices);CHKERRQ(ierr); 486 ierr = PetscFree(verticesmask);CHKERRQ(ierr); 487 ierr = PetscFree(tmp);CHKERRQ(ierr); 488 ierr = VecScatterDestroy(&user.scatter);CHKERRQ(ierr); 489 ierr = ISLocalToGlobalMappingDestroy(&isl2g);CHKERRQ(ierr); 490 ierr = VecDestroy(&x);CHKERRQ(ierr); 491 ierr = VecDestroy(&r);CHKERRQ(ierr); 492 ierr = VecDestroy(&user.localX);CHKERRQ(ierr); 493 ierr = VecDestroy(&user.localF);CHKERRQ(ierr); 494 ierr = MatDestroy(&Jac);CHKERRQ(ierr); ierr = SNESDestroy(&snes);CHKERRQ(ierr); 495 /*ierr = PetscDrawDestroy(draw);CHKERRQ(ierr);*/ 496 if (fd_jacobian_coloring) { 497 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 498 } 499 ierr = PetscFinalize(); 500 return ierr; 501 } 502 /* -------------------- Form initial approximation ----------------- */ 503 504 /* 505 FormInitialGuess - Forms initial approximation. 506 507 Input Parameters: 508 user - user-defined application context 509 X - vector 510 511 Output Parameter: 512 X - vector 513 */ 514 PetscErrorCode FormInitialGuess(AppCtx *user,Vec X) 515 { 516 PetscInt i,Nvlocal,ierr; 517 PetscInt *gloInd; 518 PetscScalar *x; 519 #if defined(UNUSED_VARIABLES) 520 PetscReal temp1,temp,hx,hy,hxdhy,hydhx,sc; 521 PetscInt Neglobal,Nvglobal,j,row; 522 PetscReal alpha,lambda; 523 524 Nvglobal = user->Nvglobal; 525 Neglobal = user->Neglobal; 526 lambda = user->non_lin_param; 527 alpha = user->lin_param; 528 #endif 529 530 Nvlocal = user->Nvlocal; 531 gloInd = user->gloInd; 532 533 /* 534 Get a pointer to vector data. 535 - For default PETSc vectors, VecGetArray() returns a pointer to 536 the data array. Otherwise, the routine is implementation dependent. 537 - You MUST call VecRestoreArray() when you no longer need access to 538 the array. 539 */ 540 ierr = VecGetArray(X,&x);CHKERRQ(ierr); 541 542 /* 543 Compute initial guess over the locally owned part of the grid 544 */ 545 for (i=0; i < Nvlocal; i++) x[i] = (PetscReal)gloInd[i]; 546 547 /* 548 Restore vector 549 */ 550 ierr = VecRestoreArray(X,&x);CHKERRQ(ierr); 551 return 0; 552 } 553 /* -------------------- Evaluate Function F(x) --------------------- */ 554 /* 555 FormFunction - Evaluates nonlinear function, F(x). 556 557 Input Parameters: 558 . snes - the SNES context 559 . X - input vector 560 . ptr - optional user-defined context, as set by SNESSetFunction() 561 562 Output Parameter: 563 . F - function vector 564 */ 565 PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr) 566 { 567 AppCtx *user = (AppCtx*)ptr; 568 PetscErrorCode ierr; 569 PetscInt i,j,Nvlocal; 570 PetscReal alpha,lambda; 571 PetscScalar *x,*f; 572 VecScatter scatter; 573 Vec localX = user->localX; 574 #if defined(UNUSED_VARIABLES) 575 PetscScalar ut,ub,ul,ur,u,*g,sc,uyy,uxx; 576 PetscReal hx,hy,hxdhy,hydhx; 577 PetscReal two = 2.0,one = 1.0; 578 PetscInt Nvglobal,Neglobal,row; 579 PetscInt *gloInd; 580 581 Nvglobal = user->Nvglobal; 582 Neglobal = user->Neglobal; 583 gloInd = user->gloInd; 584 #endif 585 586 Nvlocal = user->Nvlocal; 587 lambda = user->non_lin_param; 588 alpha = user->lin_param; 589 scatter = user->scatter; 590 591 /* 592 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 593 described in the beginning of this code 594 595 First scatter the distributed vector X into local vector localX (that includes 596 values for ghost nodes. If we wish,we can put some other work between 597 VecScatterBegin() and VecScatterEnd() to overlap the communication with 598 computation. 599 */ 600 ierr = VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 601 ierr = VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 602 603 /* 604 Get pointers to vector data 605 */ 606 ierr = VecGetArray(localX,&x);CHKERRQ(ierr); 607 ierr = VecGetArray(F,&f);CHKERRQ(ierr); 608 609 /* 610 Now compute the f(x). As mentioned earlier, the computed Laplacian is just an 611 approximate one chosen for illustrative purpose only. Another point to notice 612 is that this is a local (completly parallel) calculation. In practical application 613 codes, function calculation time is a dominat portion of the overall execution time. 614 */ 615 for (i=0; i < Nvlocal; i++) { 616 f[i] = (user->itot[i] - alpha)*x[i] - lambda*x[i]*x[i]; 617 for (j = 0; j < user->itot[i]; j++) f[i] -= x[user->AdjM[i][j]]; 618 } 619 620 /* 621 Restore vectors 622 */ 623 ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr); 624 ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 625 /*ierr = VecView(F,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr);*/ 626 627 return 0; 628 } 629 630 /* -------------------- Evaluate Jacobian F'(x) -------------------- */ 631 /* 632 FormJacobian - Evaluates Jacobian matrix. 633 634 Input Parameters: 635 . snes - the SNES context 636 . X - input vector 637 . ptr - optional user-defined context, as set by SNESSetJacobian() 638 639 Output Parameters: 640 . A - Jacobian matrix 641 . B - optionally different preconditioning matrix 642 . flag - flag indicating matrix structure 643 644 */ 645 PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr) 646 { 647 AppCtx *user = (AppCtx*)ptr; 648 PetscInt i,j,Nvlocal,col[50],ierr; 649 PetscScalar alpha,lambda,value[50]; 650 Vec localX = user->localX; 651 VecScatter scatter; 652 PetscScalar *x; 653 #if defined(UNUSED_VARIABLES) 654 PetscScalar two = 2.0,one = 1.0; 655 PetscInt row,Nvglobal,Neglobal; 656 PetscInt *gloInd; 657 658 Nvglobal = user->Nvglobal; 659 Neglobal = user->Neglobal; 660 gloInd = user->gloInd; 661 #endif 662 663 /*printf("Entering into FormJacobian \n");*/ 664 Nvlocal = user->Nvlocal; 665 lambda = user->non_lin_param; 666 alpha = user->lin_param; 667 scatter = user->scatter; 668 669 /* 670 PDE : L(u) + lambda*u*u +alpha*u = 0 where L(u) is the approximate Laplacian as 671 described in the beginning of this code 672 673 First scatter the distributed vector X into local vector localX (that includes 674 values for ghost nodes. If we wish, we can put some other work between 675 VecScatterBegin() and VecScatterEnd() to overlap the communication with 676 computation. 677 */ 678 ierr = VecScatterBegin(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 679 ierr = VecScatterEnd(scatter,X,localX,INSERT_VALUES,SCATTER_FORWARD);CHKERRQ(ierr); 680 681 /* 682 Get pointer to vector data 683 */ 684 ierr = VecGetArray(localX,&x);CHKERRQ(ierr); 685 686 for (i=0; i < Nvlocal; i++) { 687 col[0] = i; 688 value[0] = user->itot[i] - 2.0*lambda*x[i] - alpha; 689 for (j = 0; j < user->itot[i]; j++) { 690 col[j+1] = user->AdjM[i][j]; 691 value[j+1] = -1.0; 692 } 693 694 /* 695 Set the matrix values in the local ordering. Note that in order to use this 696 feature we must call the routine MatSetLocalToGlobalMapping() after the 697 matrix has been created. 698 */ 699 ierr = MatSetValuesLocal(jac,1,&i,1+user->itot[i],col,value,INSERT_VALUES);CHKERRQ(ierr); 700 } 701 702 /* 703 Assemble matrix, using the 2-step process: 704 MatAssemblyBegin(), MatAssemblyEnd(). 705 Between these two calls, the pointer to vector data has been restored to 706 demonstrate the use of overlapping communicationn with computation. 707 */ 708 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 709 ierr = VecRestoreArray(localX,&x);CHKERRQ(ierr); 710 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 711 712 /* 713 Tell the matrix we will never add a new nonzero location to the 714 matrix. If we do, it will generate an error. 715 */ 716 ierr = MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 717 /* MatView(jac,PETSC_VIEWER_STDOUT_SELF); */ 718 return 0; 719 } 720 721 /*TEST 722 723 build: 724 requires: !complex 725 726 test: 727 nsize: 2 728 args: -snes_monitor_short 729 localrunfiles: options.inf adj.in 730 731 test: 732 suffix: 2 733 nsize: 2 734 args: -snes_monitor_short -fd_jacobian_coloring 735 localrunfiles: options.inf adj.in 736 737 TEST*/ 738