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