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