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