1 2 static char help[] ="Nonlinear Radiative Transport PDE with multigrid in 2d.\n\ 3 Uses 2-dimensional distributed arrays.\n\ 4 A 2-dim simplified Radiative Transport test problem is used, with analytic Jacobian. \n\ 5 \n\ 6 Solves the linear systems via multilevel methods \n\ 7 \n\ 8 The command line\n\ 9 options are:\n\ 10 -tleft <tl>, where <tl> indicates the left Diriclet BC \n\ 11 -tright <tr>, where <tr> indicates the right Diriclet BC \n\ 12 -beta <beta>, where <beta> indicates the exponent in T \n\n"; 13 14 /* 15 16 This example models the partial differential equation 17 18 - Div(alpha* T^beta (GRAD T)) = 0. 19 20 where beta = 2.5 and alpha = 1.0 21 22 BC: T_left = 1.0, T_right = 0.1, dT/dn_top = dTdn_bottom = 0. 23 24 in the unit square, which is uniformly discretized in each of x and 25 y in this simple encoding. The degrees of freedom are cell centered. 26 27 A finite volume approximation with the usual 5-point stencil 28 is used to discretize the boundary value problem to obtain a 29 nonlinear system of equations. 30 31 This code was contributed by David Keyes 32 33 */ 34 35 #include <petscsnes.h> 36 #include <petscdm.h> 37 #include <petscdmda.h> 38 39 /* User-defined application context */ 40 41 typedef struct { 42 PetscReal tleft,tright; /* Dirichlet boundary conditions */ 43 PetscReal beta,bm1,coef; /* nonlinear diffusivity parameterizations */ 44 } AppCtx; 45 46 #define POWFLOP 5 /* assume a pow() takes five flops */ 47 48 extern PetscErrorCode FormInitialGuess(SNES,Vec,void*); 49 extern PetscErrorCode FormFunction(SNES,Vec,Vec,void*); 50 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 51 52 int main(int argc,char **argv) 53 { 54 SNES snes; 55 AppCtx user; 56 PetscInt its,lits; 57 PetscReal litspit; 58 DM da; 59 60 PetscCall(PetscInitialize(&argc,&argv,NULL,help)); 61 62 /* set problem parameters */ 63 user.tleft = 1.0; 64 user.tright = 0.1; 65 user.beta = 2.5; 66 PetscCall(PetscOptionsGetReal(NULL,NULL,"-tleft",&user.tleft,NULL)); 67 PetscCall(PetscOptionsGetReal(NULL,NULL,"-tright",&user.tright,NULL)); 68 PetscCall(PetscOptionsGetReal(NULL,NULL,"-beta",&user.beta,NULL)); 69 user.bm1 = user.beta - 1.0; 70 user.coef = user.beta/2.0; 71 72 /* 73 Create the multilevel DM data structure 74 */ 75 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes)); 76 77 /* 78 Set the DMDA (grid structure) for the grids. 79 */ 80 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_STAR,5,5,PETSC_DECIDE,PETSC_DECIDE,1,1,0,0,&da)); 81 PetscCall(DMSetFromOptions(da)); 82 PetscCall(DMSetUp(da)); 83 PetscCall(DMSetApplicationContext(da,&user)); 84 PetscCall(SNESSetDM(snes,(DM)da)); 85 86 /* 87 Create the nonlinear solver, and tell it the functions to use 88 */ 89 PetscCall(SNESSetFunction(snes,NULL,FormFunction,&user)); 90 PetscCall(SNESSetJacobian(snes,NULL,NULL,FormJacobian,&user)); 91 PetscCall(SNESSetFromOptions(snes)); 92 PetscCall(SNESSetComputeInitialGuess(snes,FormInitialGuess,NULL)); 93 94 PetscCall(SNESSolve(snes,NULL,NULL)); 95 PetscCall(SNESGetIterationNumber(snes,&its)); 96 PetscCall(SNESGetLinearSolveIterations(snes,&lits)); 97 litspit = ((PetscReal)lits)/((PetscReal)its); 98 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Number of SNES iterations = %D\n",its)); 99 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Number of Linear iterations = %D\n",lits)); 100 PetscCall(PetscPrintf(PETSC_COMM_WORLD,"Average Linear its / SNES = %e\n",(double)litspit)); 101 102 PetscCall(DMDestroy(&da)); 103 PetscCall(SNESDestroy(&snes)); 104 PetscCall(PetscFinalize()); 105 return 0; 106 } 107 /* -------------------- Form initial approximation ----------------- */ 108 PetscErrorCode FormInitialGuess(SNES snes,Vec X,void *ctx) 109 { 110 AppCtx *user; 111 PetscInt i,j,xs,ys,xm,ym; 112 PetscReal tleft; 113 PetscScalar **x; 114 DM da; 115 116 PetscFunctionBeginUser; 117 PetscCall(SNESGetDM(snes,&da)); 118 PetscCall(DMGetApplicationContext(da,&user)); 119 tleft = user->tleft; 120 /* Get ghost points */ 121 PetscCall(DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0)); 122 PetscCall(DMDAVecGetArray(da,X,&x)); 123 124 /* Compute initial guess */ 125 for (j=ys; j<ys+ym; j++) { 126 for (i=xs; i<xs+xm; i++) { 127 x[j][i] = tleft; 128 } 129 } 130 PetscCall(DMDAVecRestoreArray(da,X,&x)); 131 PetscFunctionReturn(0); 132 } 133 /* -------------------- Evaluate Function F(x) --------------------- */ 134 PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *ptr) 135 { 136 AppCtx *user = (AppCtx*)ptr; 137 PetscInt i,j,mx,my,xs,ys,xm,ym; 138 PetscScalar zero = 0.0,one = 1.0; 139 PetscScalar hx,hy,hxdhy,hydhx; 140 PetscScalar t0,tn,ts,te,tw,an,as,ae,aw,dn,ds,de,dw,fn = 0.0,fs = 0.0,fe =0.0,fw = 0.0; 141 PetscScalar tleft,tright,beta; 142 PetscScalar **x,**f; 143 Vec localX; 144 DM da; 145 146 PetscFunctionBeginUser; 147 PetscCall(SNESGetDM(snes,&da)); 148 PetscCall(DMGetLocalVector(da,&localX)); 149 PetscCall(DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0)); 150 hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); 151 hxdhy = hx/hy; hydhx = hy/hx; 152 tleft = user->tleft; tright = user->tright; 153 beta = user->beta; 154 155 /* Get ghost points */ 156 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 157 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 158 PetscCall(DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0)); 159 PetscCall(DMDAVecGetArray(da,localX,&x)); 160 PetscCall(DMDAVecGetArray(da,F,&f)); 161 162 /* Evaluate function */ 163 for (j=ys; j<ys+ym; j++) { 164 for (i=xs; i<xs+xm; i++) { 165 t0 = x[j][i]; 166 167 if (i > 0 && i < mx-1 && j > 0 && j < my-1) { 168 169 /* general interior volume */ 170 171 tw = x[j][i-1]; 172 aw = 0.5*(t0 + tw); 173 dw = PetscPowScalar(aw,beta); 174 fw = dw*(t0 - tw); 175 176 te = x[j][i+1]; 177 ae = 0.5*(t0 + te); 178 de = PetscPowScalar(ae,beta); 179 fe = de*(te - t0); 180 181 ts = x[j-1][i]; 182 as = 0.5*(t0 + ts); 183 ds = PetscPowScalar(as,beta); 184 fs = ds*(t0 - ts); 185 186 tn = x[j+1][i]; 187 an = 0.5*(t0 + tn); 188 dn = PetscPowScalar(an,beta); 189 fn = dn*(tn - t0); 190 191 } else if (i == 0) { 192 193 /* left-hand boundary */ 194 tw = tleft; 195 aw = 0.5*(t0 + tw); 196 dw = PetscPowScalar(aw,beta); 197 fw = dw*(t0 - tw); 198 199 te = x[j][i+1]; 200 ae = 0.5*(t0 + te); 201 de = PetscPowScalar(ae,beta); 202 fe = de*(te - t0); 203 204 if (j > 0) { 205 ts = x[j-1][i]; 206 as = 0.5*(t0 + ts); 207 ds = PetscPowScalar(as,beta); 208 fs = ds*(t0 - ts); 209 } else fs = zero; 210 211 if (j < my-1) { 212 tn = x[j+1][i]; 213 an = 0.5*(t0 + tn); 214 dn = PetscPowScalar(an,beta); 215 fn = dn*(tn - t0); 216 } else fn = zero; 217 218 } else if (i == mx-1) { 219 220 /* right-hand boundary */ 221 tw = x[j][i-1]; 222 aw = 0.5*(t0 + tw); 223 dw = PetscPowScalar(aw,beta); 224 fw = dw*(t0 - tw); 225 226 te = tright; 227 ae = 0.5*(t0 + te); 228 de = PetscPowScalar(ae,beta); 229 fe = de*(te - t0); 230 231 if (j > 0) { 232 ts = x[j-1][i]; 233 as = 0.5*(t0 + ts); 234 ds = PetscPowScalar(as,beta); 235 fs = ds*(t0 - ts); 236 } else fs = zero; 237 238 if (j < my-1) { 239 tn = x[j+1][i]; 240 an = 0.5*(t0 + tn); 241 dn = PetscPowScalar(an,beta); 242 fn = dn*(tn - t0); 243 } else fn = zero; 244 245 } else if (j == 0) { 246 247 /* bottom boundary,and i <> 0 or mx-1 */ 248 tw = x[j][i-1]; 249 aw = 0.5*(t0 + tw); 250 dw = PetscPowScalar(aw,beta); 251 fw = dw*(t0 - tw); 252 253 te = x[j][i+1]; 254 ae = 0.5*(t0 + te); 255 de = PetscPowScalar(ae,beta); 256 fe = de*(te - t0); 257 258 fs = zero; 259 260 tn = x[j+1][i]; 261 an = 0.5*(t0 + tn); 262 dn = PetscPowScalar(an,beta); 263 fn = dn*(tn - t0); 264 265 } else if (j == my-1) { 266 267 /* top boundary,and i <> 0 or mx-1 */ 268 tw = x[j][i-1]; 269 aw = 0.5*(t0 + tw); 270 dw = PetscPowScalar(aw,beta); 271 fw = dw*(t0 - tw); 272 273 te = x[j][i+1]; 274 ae = 0.5*(t0 + te); 275 de = PetscPowScalar(ae,beta); 276 fe = de*(te - t0); 277 278 ts = x[j-1][i]; 279 as = 0.5*(t0 + ts); 280 ds = PetscPowScalar(as,beta); 281 fs = ds*(t0 - ts); 282 283 fn = zero; 284 285 } 286 287 f[j][i] = -hydhx*(fe-fw) - hxdhy*(fn-fs); 288 289 } 290 } 291 PetscCall(DMDAVecRestoreArray(da,localX,&x)); 292 PetscCall(DMDAVecRestoreArray(da,F,&f)); 293 PetscCall(DMRestoreLocalVector(da,&localX)); 294 PetscCall(PetscLogFlops((22.0 + 4.0*POWFLOP)*ym*xm)); 295 PetscFunctionReturn(0); 296 } 297 /* -------------------- Evaluate Jacobian F(x) --------------------- */ 298 PetscErrorCode FormJacobian(SNES snes,Vec X,Mat jac,Mat B,void *ptr) 299 { 300 AppCtx *user = (AppCtx*)ptr; 301 PetscInt i,j,mx,my,xs,ys,xm,ym; 302 PetscScalar one = 1.0,hx,hy,hxdhy,hydhx,t0,tn,ts,te,tw; 303 PetscScalar dn,ds,de,dw,an,as,ae,aw,bn,bs,be,bw,gn,gs,ge,gw; 304 PetscScalar tleft,tright,beta,bm1,coef; 305 PetscScalar v[5],**x; 306 Vec localX; 307 MatStencil col[5],row; 308 DM da; 309 310 PetscFunctionBeginUser; 311 PetscCall(SNESGetDM(snes,&da)); 312 PetscCall(DMGetLocalVector(da,&localX)); 313 PetscCall(DMDAGetInfo(da,NULL,&mx,&my,0,0,0,0,0,0,0,0,0,0)); 314 hx = one/(PetscReal)(mx-1); hy = one/(PetscReal)(my-1); 315 hxdhy = hx/hy; hydhx = hy/hx; 316 tleft = user->tleft; tright = user->tright; 317 beta = user->beta; bm1 = user->bm1; coef = user->coef; 318 319 /* Get ghost points */ 320 PetscCall(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 321 PetscCall(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 322 PetscCall(DMDAGetCorners(da,&xs,&ys,0,&xm,&ym,0)); 323 PetscCall(DMDAVecGetArray(da,localX,&x)); 324 325 /* Evaluate Jacobian of function */ 326 for (j=ys; j<ys+ym; j++) { 327 for (i=xs; i<xs+xm; i++) { 328 t0 = x[j][i]; 329 330 if (i > 0 && i < mx-1 && j > 0 && j < my-1) { 331 332 /* general interior volume */ 333 334 tw = x[j][i-1]; 335 aw = 0.5*(t0 + tw); 336 bw = PetscPowScalar(aw,bm1); 337 /* dw = bw * aw */ 338 dw = PetscPowScalar(aw,beta); 339 gw = coef*bw*(t0 - tw); 340 341 te = x[j][i+1]; 342 ae = 0.5*(t0 + te); 343 be = PetscPowScalar(ae,bm1); 344 /* de = be * ae; */ 345 de = PetscPowScalar(ae,beta); 346 ge = coef*be*(te - t0); 347 348 ts = x[j-1][i]; 349 as = 0.5*(t0 + ts); 350 bs = PetscPowScalar(as,bm1); 351 /* ds = bs * as; */ 352 ds = PetscPowScalar(as,beta); 353 gs = coef*bs*(t0 - ts); 354 355 tn = x[j+1][i]; 356 an = 0.5*(t0 + tn); 357 bn = PetscPowScalar(an,bm1); 358 /* dn = bn * an; */ 359 dn = PetscPowScalar(an,beta); 360 gn = coef*bn*(tn - t0); 361 362 v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 363 v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 364 v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 365 v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1; 366 v[4] = -hxdhy*(dn + gn); col[4].j = j+1; col[4].i = i; 367 PetscCall(MatSetValuesStencil(B,1,&row,5,col,v,INSERT_VALUES)); 368 369 } else if (i == 0) { 370 371 /* left-hand boundary */ 372 tw = tleft; 373 aw = 0.5*(t0 + tw); 374 bw = PetscPowScalar(aw,bm1); 375 /* dw = bw * aw */ 376 dw = PetscPowScalar(aw,beta); 377 gw = coef*bw*(t0 - tw); 378 379 te = x[j][i + 1]; 380 ae = 0.5*(t0 + te); 381 be = PetscPowScalar(ae,bm1); 382 /* de = be * ae; */ 383 de = PetscPowScalar(ae,beta); 384 ge = coef*be*(te - t0); 385 386 /* left-hand bottom boundary */ 387 if (j == 0) { 388 389 tn = x[j+1][i]; 390 an = 0.5*(t0 + tn); 391 bn = PetscPowScalar(an,bm1); 392 /* dn = bn * an; */ 393 dn = PetscPowScalar(an,beta); 394 gn = coef*bn*(tn - t0); 395 396 v[0] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[0].j = row.j = j; col[0].i = row.i = i; 397 v[1] = -hydhx*(de + ge); col[1].j = j; col[1].i = i+1; 398 v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i; 399 PetscCall(MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES)); 400 401 /* left-hand interior boundary */ 402 } else if (j < my-1) { 403 404 ts = x[j-1][i]; 405 as = 0.5*(t0 + ts); 406 bs = PetscPowScalar(as,bm1); 407 /* ds = bs * as; */ 408 ds = PetscPowScalar(as,beta); 409 gs = coef*bs*(t0 - ts); 410 411 tn = x[j+1][i]; 412 an = 0.5*(t0 + tn); 413 bn = PetscPowScalar(an,bm1); 414 /* dn = bn * an; */ 415 dn = PetscPowScalar(an,beta); 416 gn = coef*bn*(tn - t0); 417 418 v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 419 v[1] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 420 v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 421 v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 422 PetscCall(MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES)); 423 /* left-hand top boundary */ 424 } else { 425 426 ts = x[j-1][i]; 427 as = 0.5*(t0 + ts); 428 bs = PetscPowScalar(as,bm1); 429 /* ds = bs * as; */ 430 ds = PetscPowScalar(as,beta); 431 gs = coef*bs*(t0 - ts); 432 433 v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 434 v[1] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 435 v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 436 PetscCall(MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES)); 437 } 438 439 } else if (i == mx-1) { 440 441 /* right-hand boundary */ 442 tw = x[j][i-1]; 443 aw = 0.5*(t0 + tw); 444 bw = PetscPowScalar(aw,bm1); 445 /* dw = bw * aw */ 446 dw = PetscPowScalar(aw,beta); 447 gw = coef*bw*(t0 - tw); 448 449 te = tright; 450 ae = 0.5*(t0 + te); 451 be = PetscPowScalar(ae,bm1); 452 /* de = be * ae; */ 453 de = PetscPowScalar(ae,beta); 454 ge = coef*be*(te - t0); 455 456 /* right-hand bottom boundary */ 457 if (j == 0) { 458 459 tn = x[j+1][i]; 460 an = 0.5*(t0 + tn); 461 bn = PetscPowScalar(an,bm1); 462 /* dn = bn * an; */ 463 dn = PetscPowScalar(an,beta); 464 gn = coef*bn*(tn - t0); 465 466 v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1; 467 v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 468 v[2] = -hxdhy*(dn + gn); col[2].j = j+1; col[2].i = i; 469 PetscCall(MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES)); 470 471 /* right-hand interior boundary */ 472 } else if (j < my-1) { 473 474 ts = x[j-1][i]; 475 as = 0.5*(t0 + ts); 476 bs = PetscPowScalar(as,bm1); 477 /* ds = bs * as; */ 478 ds = PetscPowScalar(as,beta); 479 gs = coef*bs*(t0 - ts); 480 481 tn = x[j+1][i]; 482 an = 0.5*(t0 + tn); 483 bn = PetscPowScalar(an,bm1); 484 /* dn = bn * an; */ 485 dn = PetscPowScalar(an,beta); 486 gn = coef*bn*(tn - t0); 487 488 v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 489 v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 490 v[2] = hxdhy*(ds + dn + gs - gn) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 491 v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 492 PetscCall(MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES)); 493 /* right-hand top boundary */ 494 } else { 495 496 ts = x[j-1][i]; 497 as = 0.5*(t0 + ts); 498 bs = PetscPowScalar(as,bm1); 499 /* ds = bs * as; */ 500 ds = PetscPowScalar(as,beta); 501 gs = coef*bs*(t0 - ts); 502 503 v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 504 v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 505 v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 506 PetscCall(MatSetValuesStencil(B,1,&row,3,col,v,INSERT_VALUES)); 507 } 508 509 /* bottom boundary,and i <> 0 or mx-1 */ 510 } else if (j == 0) { 511 512 tw = x[j][i-1]; 513 aw = 0.5*(t0 + tw); 514 bw = PetscPowScalar(aw,bm1); 515 /* dw = bw * aw */ 516 dw = PetscPowScalar(aw,beta); 517 gw = coef*bw*(t0 - tw); 518 519 te = x[j][i+1]; 520 ae = 0.5*(t0 + te); 521 be = PetscPowScalar(ae,bm1); 522 /* de = be * ae; */ 523 de = PetscPowScalar(ae,beta); 524 ge = coef*be*(te - t0); 525 526 tn = x[j+1][i]; 527 an = 0.5*(t0 + tn); 528 bn = PetscPowScalar(an,bm1); 529 /* dn = bn * an; */ 530 dn = PetscPowScalar(an,beta); 531 gn = coef*bn*(tn - t0); 532 533 v[0] = -hydhx*(dw - gw); col[0].j = j; col[0].i = i-1; 534 v[1] = hxdhy*(dn - gn) + hydhx*(dw + de + gw - ge); col[1].j = row.j = j; col[1].i = row.i = i; 535 v[2] = -hydhx*(de + ge); col[2].j = j; col[2].i = i+1; 536 v[3] = -hxdhy*(dn + gn); col[3].j = j+1; col[3].i = i; 537 PetscCall(MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES)); 538 539 /* top boundary,and i <> 0 or mx-1 */ 540 } else if (j == my-1) { 541 542 tw = x[j][i-1]; 543 aw = 0.5*(t0 + tw); 544 bw = PetscPowScalar(aw,bm1); 545 /* dw = bw * aw */ 546 dw = PetscPowScalar(aw,beta); 547 gw = coef*bw*(t0 - tw); 548 549 te = x[j][i+1]; 550 ae = 0.5*(t0 + te); 551 be = PetscPowScalar(ae,bm1); 552 /* de = be * ae; */ 553 de = PetscPowScalar(ae,beta); 554 ge = coef*be*(te - t0); 555 556 ts = x[j-1][i]; 557 as = 0.5*(t0 + ts); 558 bs = PetscPowScalar(as,bm1); 559 /* ds = bs * as; */ 560 ds = PetscPowScalar(as,beta); 561 gs = coef*bs*(t0 - ts); 562 563 v[0] = -hxdhy*(ds - gs); col[0].j = j-1; col[0].i = i; 564 v[1] = -hydhx*(dw - gw); col[1].j = j; col[1].i = i-1; 565 v[2] = hxdhy*(ds + gs) + hydhx*(dw + de + gw - ge); col[2].j = row.j = j; col[2].i = row.i = i; 566 v[3] = -hydhx*(de + ge); col[3].j = j; col[3].i = i+1; 567 PetscCall(MatSetValuesStencil(B,1,&row,4,col,v,INSERT_VALUES)); 568 569 } 570 } 571 } 572 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 573 PetscCall(DMDAVecRestoreArray(da,localX,&x)); 574 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 575 PetscCall(DMRestoreLocalVector(da,&localX)); 576 if (jac != B) { 577 PetscCall(MatAssemblyBegin(jac,MAT_FINAL_ASSEMBLY)); 578 PetscCall(MatAssemblyEnd(jac,MAT_FINAL_ASSEMBLY)); 579 } 580 581 PetscCall(PetscLogFlops((41.0 + 8.0*POWFLOP)*xm*ym)); 582 PetscFunctionReturn(0); 583 } 584 585 /*TEST 586 587 test: 588 suffix: 1 589 args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view 590 requires: !single 591 592 test: 593 suffix: 2 594 args: -pc_type mg -ksp_type fgmres -da_refine 2 -pc_mg_galerkin pmat -snes_view -snes_type newtontrdc 595 requires: !single 596 597 TEST*/ 598