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