1 2 static char help[] = "Bratu nonlinear PDE in 3d.\n\ 3 We solve the Bratu (SFI - solid fuel ignition) problem in a 3D rectangular\n\ 4 domain, using distributed arrays (DMDAs) to partition the parallel grid.\n\ 5 The command line options include:\n\ 6 -par <parameter>, where <parameter> indicates the problem's nonlinearity\n\ 7 problem SFI: <parameter> = Bratu parameter (0 <= par <= 6.81)\n\n"; 8 9 /*T 10 Concepts: SNES^parallel Bratu example 11 Concepts: DMDA^using distributed arrays; 12 Processors: n 13 T*/ 14 15 16 17 /* ------------------------------------------------------------------------ 18 19 Solid Fuel Ignition (SFI) problem. This problem is modeled by 20 the partial differential equation 21 22 -Laplacian u - lambda*exp(u) = 0, 0 < x,y < 1, 23 24 with boundary conditions 25 26 u = 0 for x = 0, x = 1, y = 0, y = 1, z = 0, z = 1 27 28 A finite difference approximation with the usual 7-point stencil 29 is used to discretize the boundary value problem to obtain a nonlinear 30 system of equations. 31 32 33 ------------------------------------------------------------------------- */ 34 35 /* 36 Include "petscdmda.h" so that we can use distributed arrays (DMDAs). 37 Include "petscsnes.h" so that we can use SNES solvers. Note that this 38 file automatically includes: 39 petscsys.h - base PETSc routines petscvec.h - vectors 40 petscmat.h - matrices 41 petscis.h - index sets petscksp.h - Krylov subspace methods 42 petscviewer.h - viewers petscpc.h - preconditioners 43 petscksp.h - linear solvers 44 */ 45 #include <petscdm.h> 46 #include <petscdmda.h> 47 #include <petscsnes.h> 48 49 50 /* 51 User-defined application context - contains data needed by the 52 application-provided call-back routines, FormJacobian() and 53 FormFunction(). 54 */ 55 typedef struct { 56 PetscReal param; /* test problem parameter */ 57 DM da; /* distributed array data structure */ 58 } AppCtx; 59 60 /* 61 User-defined routines 62 */ 63 extern PetscErrorCode FormFunctionLocal(SNES,Vec,Vec,void*); 64 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 65 extern PetscErrorCode FormInitialGuess(AppCtx*,Vec); 66 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 67 68 int main(int argc,char **argv) 69 { 70 SNES snes; /* nonlinear solver */ 71 Vec x,r; /* solution, residual vectors */ 72 Mat J = NULL; /* Jacobian matrix */ 73 AppCtx user; /* user-defined work context */ 74 PetscInt its; /* iterations for convergence */ 75 MatFDColoring matfdcoloring = NULL; 76 PetscBool matrix_free = PETSC_FALSE,coloring = PETSC_FALSE, coloring_ds = PETSC_FALSE,local_coloring = PETSC_FALSE; 77 PetscErrorCode ierr; 78 PetscReal bratu_lambda_max = 6.81,bratu_lambda_min = 0.,fnorm; 79 80 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 81 Initialize program 82 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 83 84 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 85 86 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 87 Initialize problem parameters 88 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 89 user.param = 6.0; 90 ierr = PetscOptionsGetReal(NULL,NULL,"-par",&user.param,NULL);CHKERRQ(ierr); 91 if (user.param >= bratu_lambda_max || user.param <= bratu_lambda_min) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Lambda is out of range"); 92 93 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 94 Create nonlinear solver context 95 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 96 ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr); 97 98 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 99 Create distributed array (DMDA) to manage parallel grid and vectors 100 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 101 ierr = DMDACreate3d(PETSC_COMM_WORLD,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,4,4,4,PETSC_DECIDE,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,NULL,&user.da);CHKERRQ(ierr); 102 ierr = DMSetFromOptions(user.da);CHKERRQ(ierr); 103 ierr = DMSetUp(user.da);CHKERRQ(ierr); 104 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 105 Extract global vectors from DMDA; then duplicate for remaining 106 vectors that are the same types 107 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 108 ierr = DMCreateGlobalVector(user.da,&x);CHKERRQ(ierr); 109 ierr = VecDuplicate(x,&r);CHKERRQ(ierr); 110 111 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 112 Set function evaluation routine and vector 113 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 114 ierr = SNESSetFunction(snes,r,FormFunction,(void*)&user);CHKERRQ(ierr); 115 116 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 117 Create matrix data structure; set Jacobian evaluation routine 118 119 Set Jacobian matrix data structure and default Jacobian evaluation 120 routine. User can override with: 121 -snes_mf : matrix-free Newton-Krylov method with no preconditioning 122 (unless user explicitly sets preconditioner) 123 -snes_mf_operator : form preconditioning matrix as set by the user, 124 but use matrix-free approx for Jacobian-vector 125 products within Newton-Krylov method 126 -fdcoloring : using finite differences with coloring to compute the Jacobian 127 128 Note one can use -matfd_coloring wp or ds the only reason for the -fdcoloring_ds option 129 below is to test the call to MatFDColoringSetType(). 130 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 131 ierr = PetscOptionsGetBool(NULL,NULL,"-snes_mf",&matrix_free,NULL);CHKERRQ(ierr); 132 ierr = PetscOptionsGetBool(NULL,NULL,"-fdcoloring",&coloring,NULL);CHKERRQ(ierr); 133 ierr = PetscOptionsGetBool(NULL,NULL,"-fdcoloring_ds",&coloring_ds,NULL);CHKERRQ(ierr); 134 ierr = PetscOptionsGetBool(NULL,NULL,"-fdcoloring_local",&local_coloring,NULL);CHKERRQ(ierr); 135 if (!matrix_free) { 136 ierr = DMSetMatType(user.da,MATAIJ);CHKERRQ(ierr); 137 ierr = DMCreateMatrix(user.da,&J);CHKERRQ(ierr); 138 if (coloring) { 139 ISColoring iscoloring; 140 if (!local_coloring) { 141 ierr = DMCreateColoring(user.da,IS_COLORING_GLOBAL,&iscoloring);CHKERRQ(ierr); 142 ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr); 143 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunction,&user);CHKERRQ(ierr); 144 } else { 145 ierr = DMCreateColoring(user.da,IS_COLORING_LOCAL,&iscoloring);CHKERRQ(ierr); 146 ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr); 147 ierr = MatFDColoringUseDM(J,matfdcoloring);CHKERRQ(ierr); 148 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))FormFunctionLocal,&user);CHKERRQ(ierr); 149 } 150 if (coloring_ds) { 151 ierr = MatFDColoringSetType(matfdcoloring,MATMFFD_DS);CHKERRQ(ierr); 152 } 153 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr); 154 ierr = MatFDColoringSetUp(J,iscoloring,matfdcoloring);CHKERRQ(ierr); 155 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr); 156 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr); 157 } else { 158 ierr = SNESSetJacobian(snes,J,J,FormJacobian,&user);CHKERRQ(ierr); 159 } 160 } 161 162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 163 Customize nonlinear solver; set runtime options 164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 165 ierr = SNESSetDM(snes,user.da);CHKERRQ(ierr); 166 ierr = SNESSetFromOptions(snes);CHKERRQ(ierr); 167 168 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 169 Evaluate initial guess 170 Note: The user should initialize the vector, x, with the initial guess 171 for the nonlinear solver prior to calling SNESSolve(). In particular, 172 to employ an initial guess of zero, the user should explicitly set 173 this vector to zero by calling VecSet(). 174 */ 175 ierr = FormInitialGuess(&user,x);CHKERRQ(ierr); 176 177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178 Solve nonlinear system 179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 180 ierr = SNESSolve(snes,NULL,x);CHKERRQ(ierr); 181 ierr = SNESGetIterationNumber(snes,&its);CHKERRQ(ierr); 182 183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 Explicitly check norm of the residual of the solution 185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 186 ierr = FormFunction(snes,x,r,(void*)&user);CHKERRQ(ierr); 187 ierr = VecNorm(r,NORM_2,&fnorm);CHKERRQ(ierr); 188 ierr = PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D fnorm %g\n",its,(double)fnorm);CHKERRQ(ierr); 189 190 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 191 Free work space. All PETSc objects should be destroyed when they 192 are no longer needed. 193 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 194 195 ierr = MatDestroy(&J);CHKERRQ(ierr); 196 ierr = VecDestroy(&x);CHKERRQ(ierr); 197 ierr = VecDestroy(&r);CHKERRQ(ierr); 198 ierr = SNESDestroy(&snes);CHKERRQ(ierr); 199 ierr = DMDestroy(&user.da);CHKERRQ(ierr); 200 ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr); 201 ierr = PetscFinalize(); 202 return ierr; 203 } 204 /* ------------------------------------------------------------------- */ 205 /* 206 FormInitialGuess - Forms initial approximation. 207 208 Input Parameters: 209 user - user-defined application context 210 X - vector 211 212 Output Parameter: 213 X - vector 214 */ 215 PetscErrorCode FormInitialGuess(AppCtx *user,Vec X) 216 { 217 PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm; 218 PetscErrorCode ierr; 219 PetscReal lambda,temp1,hx,hy,hz,tempk,tempj; 220 PetscScalar ***x; 221 222 PetscFunctionBeginUser; 223 ierr = DMDAGetInfo(user->da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 224 225 lambda = user->param; 226 hx = 1.0/(PetscReal)(Mx-1); 227 hy = 1.0/(PetscReal)(My-1); 228 hz = 1.0/(PetscReal)(Mz-1); 229 temp1 = lambda/(lambda + 1.0); 230 231 /* 232 Get a pointer to vector data. 233 - For default PETSc vectors, VecGetArray() returns a pointer to 234 the data array. Otherwise, the routine is implementation dependent. 235 - You MUST call VecRestoreArray() when you no longer need access to 236 the array. 237 */ 238 ierr = DMDAVecGetArray(user->da,X,&x);CHKERRQ(ierr); 239 240 /* 241 Get local grid boundaries (for 3-dimensional DMDA): 242 xs, ys, zs - starting grid indices (no ghost points) 243 xm, ym, zm - widths of local grid (no ghost points) 244 245 */ 246 ierr = DMDAGetCorners(user->da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 247 248 /* 249 Compute initial guess over the locally owned part of the grid 250 */ 251 for (k=zs; k<zs+zm; k++) { 252 tempk = (PetscReal)(PetscMin(k,Mz-k-1))*hz; 253 for (j=ys; j<ys+ym; j++) { 254 tempj = PetscMin((PetscReal)(PetscMin(j,My-j-1))*hy,tempk); 255 for (i=xs; i<xs+xm; i++) { 256 if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) { 257 /* boundary conditions are all zero Dirichlet */ 258 x[k][j][i] = 0.0; 259 } else { 260 x[k][j][i] = temp1*PetscSqrtReal(PetscMin((PetscReal)(PetscMin(i,Mx-i-1))*hx,tempj)); 261 } 262 } 263 } 264 } 265 266 /* 267 Restore vector 268 */ 269 ierr = DMDAVecRestoreArray(user->da,X,&x);CHKERRQ(ierr); 270 PetscFunctionReturn(0); 271 } 272 /* ------------------------------------------------------------------- */ 273 /* 274 FormFunctionLocal - Evaluates nonlinear function, F(x) on a ghosted patch 275 276 Input Parameters: 277 . snes - the SNES context 278 . localX - input vector, this contains the ghosted region needed 279 . ptr - optional user-defined context, as set by SNESSetFunction() 280 281 Output Parameter: 282 . F - function vector, this does not contain a ghosted region 283 */ 284 PetscErrorCode FormFunctionLocal(SNES snes,Vec localX,Vec F,void *ptr) 285 { 286 AppCtx *user = (AppCtx*)ptr; 287 PetscErrorCode ierr; 288 PetscInt i,j,k,Mx,My,Mz,xs,ys,zs,xm,ym,zm; 289 PetscReal two = 2.0,lambda,hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc; 290 PetscScalar u_north,u_south,u_east,u_west,u_up,u_down,u,u_xx,u_yy,u_zz,***x,***f; 291 DM da; 292 293 PetscFunctionBeginUser; 294 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 295 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 296 297 lambda = user->param; 298 hx = 1.0/(PetscReal)(Mx-1); 299 hy = 1.0/(PetscReal)(My-1); 300 hz = 1.0/(PetscReal)(Mz-1); 301 sc = hx*hy*hz*lambda; 302 hxhzdhy = hx*hz/hy; 303 hyhzdhx = hy*hz/hx; 304 hxhydhz = hx*hy/hz; 305 306 /* 307 Get pointers to vector data 308 */ 309 ierr = DMDAVecGetArrayRead(da,localX,&x);CHKERRQ(ierr); 310 ierr = DMDAVecGetArray(da,F,&f);CHKERRQ(ierr); 311 312 /* 313 Get local grid boundaries 314 */ 315 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 316 317 /* 318 Compute function over the locally owned part of the grid 319 */ 320 for (k=zs; k<zs+zm; k++) { 321 for (j=ys; j<ys+ym; j++) { 322 for (i=xs; i<xs+xm; i++) { 323 if (i == 0 || j == 0 || k == 0 || i == Mx-1 || j == My-1 || k == Mz-1) { 324 f[k][j][i] = x[k][j][i]; 325 } else { 326 u = x[k][j][i]; 327 u_east = x[k][j][i+1]; 328 u_west = x[k][j][i-1]; 329 u_north = x[k][j+1][i]; 330 u_south = x[k][j-1][i]; 331 u_up = x[k+1][j][i]; 332 u_down = x[k-1][j][i]; 333 u_xx = (-u_east + two*u - u_west)*hyhzdhx; 334 u_yy = (-u_north + two*u - u_south)*hxhzdhy; 335 u_zz = (-u_up + two*u - u_down)*hxhydhz; 336 f[k][j][i] = u_xx + u_yy + u_zz - sc*PetscExpScalar(u); 337 } 338 } 339 } 340 } 341 342 /* 343 Restore vectors 344 */ 345 ierr = DMDAVecRestoreArrayRead(da,localX,&x);CHKERRQ(ierr); 346 ierr = DMDAVecRestoreArray(da,F,&f);CHKERRQ(ierr); 347 ierr = PetscLogFlops(11.0*ym*xm);CHKERRQ(ierr); 348 PetscFunctionReturn(0); 349 } 350 /* ------------------------------------------------------------------- */ 351 /* 352 FormFunction - Evaluates nonlinear function, F(x) on the entire domain 353 354 Input Parameters: 355 . snes - the SNES context 356 . X - input vector 357 . ptr - optional user-defined context, as set by SNESSetFunction() 358 359 Output Parameter: 360 . F - function vector 361 */ 362 PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr) 363 { 364 PetscErrorCode ierr; 365 Vec localX; 366 DM da; 367 368 PetscFunctionBeginUser; 369 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 370 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 371 372 /* 373 Scatter ghost points to local vector,using the 2-step process 374 DMGlobalToLocalBegin(),DMGlobalToLocalEnd(). 375 By placing code between these two statements, computations can be 376 done while messages are in transition. 377 */ 378 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 379 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 380 381 ierr = FormFunctionLocal(snes,localX,F,ptr);CHKERRQ(ierr); 382 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 383 PetscFunctionReturn(0); 384 } 385 /* ------------------------------------------------------------------- */ 386 /* 387 FormJacobian - Evaluates Jacobian matrix. 388 389 Input Parameters: 390 . snes - the SNES context 391 . x - input vector 392 . ptr - optional user-defined context, as set by SNESSetJacobian() 393 394 Output Parameters: 395 . A - Jacobian matrix 396 . B - optionally different preconditioning matrix 397 398 */ 399 PetscErrorCode FormJacobian(SNES snes,Vec X,Mat J,Mat jac,void *ptr) 400 { 401 AppCtx *user = (AppCtx*)ptr; /* user-defined application context */ 402 Vec localX; 403 PetscErrorCode ierr; 404 PetscInt i,j,k,Mx,My,Mz; 405 MatStencil col[7],row; 406 PetscInt xs,ys,zs,xm,ym,zm; 407 PetscScalar lambda,v[7],hx,hy,hz,hxhzdhy,hyhzdhx,hxhydhz,sc,***x; 408 DM da; 409 410 PetscFunctionBeginUser; 411 ierr = SNESGetDM(snes,&da);CHKERRQ(ierr); 412 ierr = DMGetLocalVector(da,&localX);CHKERRQ(ierr); 413 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,&Mz,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr); 414 415 lambda = user->param; 416 hx = 1.0/(PetscReal)(Mx-1); 417 hy = 1.0/(PetscReal)(My-1); 418 hz = 1.0/(PetscReal)(Mz-1); 419 sc = hx*hy*hz*lambda; 420 hxhzdhy = hx*hz/hy; 421 hyhzdhx = hy*hz/hx; 422 hxhydhz = hx*hy/hz; 423 424 /* 425 Scatter ghost points to local vector, using the 2-step process 426 DMGlobalToLocalBegin(), DMGlobalToLocalEnd(). 427 By placing code between these two statements, computations can be 428 done while messages are in transition. 429 */ 430 ierr = DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 431 ierr = DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX);CHKERRQ(ierr); 432 433 /* 434 Get pointer to vector data 435 */ 436 ierr = DMDAVecGetArrayRead(da,localX,&x);CHKERRQ(ierr); 437 438 /* 439 Get local grid boundaries 440 */ 441 ierr = DMDAGetCorners(da,&xs,&ys,&zs,&xm,&ym,&zm);CHKERRQ(ierr); 442 443 /* 444 Compute entries for the locally owned part of the Jacobian. 445 - Currently, all PETSc parallel matrix formats are partitioned by 446 contiguous chunks of rows across the processors. 447 - Each processor needs to insert only elements that it owns 448 locally (but any non-local elements will be sent to the 449 appropriate processor during matrix assembly). 450 - Here, we set all entries for a particular row at once. 451 - We can set matrix entries either using either 452 MatSetValuesLocal() or MatSetValues(), as discussed above. 453 */ 454 for (k=zs; k<zs+zm; k++) { 455 for (j=ys; j<ys+ym; j++) { 456 for (i=xs; i<xs+xm; i++) { 457 row.k = k; row.j = j; row.i = i; 458 /* boundary points */ 459 if (i == 0 || j == 0 || k == 0|| i == Mx-1 || j == My-1 || k == Mz-1) { 460 v[0] = 1.0; 461 ierr = MatSetValuesStencil(jac,1,&row,1,&row,v,INSERT_VALUES);CHKERRQ(ierr); 462 } else { 463 /* interior grid points */ 464 v[0] = -hxhydhz; col[0].k=k-1;col[0].j=j; col[0].i = i; 465 v[1] = -hxhzdhy; col[1].k=k; col[1].j=j-1;col[1].i = i; 466 v[2] = -hyhzdhx; col[2].k=k; col[2].j=j; col[2].i = i-1; 467 v[3] = 2.0*(hyhzdhx+hxhzdhy+hxhydhz)-sc*PetscExpScalar(x[k][j][i]);col[3].k=row.k;col[3].j=row.j;col[3].i = row.i; 468 v[4] = -hyhzdhx; col[4].k=k; col[4].j=j; col[4].i = i+1; 469 v[5] = -hxhzdhy; col[5].k=k; col[5].j=j+1;col[5].i = i; 470 v[6] = -hxhydhz; col[6].k=k+1;col[6].j=j; col[6].i = i; 471 ierr = MatSetValuesStencil(jac,1,&row,7,col,v,INSERT_VALUES);CHKERRQ(ierr); 472 } 473 } 474 } 475 } 476 ierr = DMDAVecRestoreArrayRead(da,localX,&x);CHKERRQ(ierr); 477 ierr = DMRestoreLocalVector(da,&localX);CHKERRQ(ierr); 478 479 /* 480 Assemble matrix, using the 2-step process: 481 MatAssemblyBegin(), MatAssemblyEnd(). 482 */ 483 ierr = MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 484 ierr = MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 485 486 /* 487 Normally since the matrix has already been assembled above; this 488 would do nothing. But in the matrix free mode -snes_mf_operator 489 this tells the "matrix-free" matrix that a new linear system solve 490 is about to be done. 491 */ 492 493 ierr = MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 494 ierr = MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 495 496 /* 497 Tell the matrix we will never add a new nonzero location to the 498 matrix. If we do, it will generate an error. 499 */ 500 ierr = MatSetOption(jac,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);CHKERRQ(ierr); 501 PetscFunctionReturn(0); 502 } 503 504 505 506 /*TEST 507 508 test: 509 nsize: 4 510 args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 511 512 test: 513 suffix: 2 514 nsize: 4 515 args: -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 516 517 test: 518 suffix: 3 519 nsize: 4 520 args: -fdcoloring -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 521 522 test: 523 suffix: 3_ds 524 nsize: 4 525 args: -fdcoloring -fdcoloring_ds -snes_monitor_short -ksp_gmres_cgs_refinement_type refine_always 526 527 test: 528 suffix: 4 529 nsize: 4 530 args: -fdcoloring_local -fdcoloring -ksp_monitor_short -da_refine 1 531 requires: !single 532 533 TEST*/ 534