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