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