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