1 2 #include <petscsnes.h> 3 #include <petscdm.h> 4 #include <petscdmda.h> 5 6 static const char help[] = "Parallel version of the minimum surface area problem in 2D using DMDA.\n\ 7 It solves a system of nonlinear equations in mixed\n\ 8 complementarity form.This example is based on a\n\ 9 problem from the MINPACK-2 test suite. Given a rectangular 2-D domain and\n\ 10 boundary values along the edges of the domain, the objective is to find the\n\ 11 surface with the minimal area that satisfies the boundary conditions.\n\ 12 This application solves this problem using complimentarity -- We are actually\n\ 13 solving the system (grad f)_i >= 0, if x_i == l_i \n\ 14 (grad f)_i = 0, if l_i < x_i < u_i \n\ 15 (grad f)_i <= 0, if x_i == u_i \n\ 16 where f is the function to be minimized. \n\ 17 \n\ 18 The command line options are:\n\ 19 -da_grid_x <nx>, where <nx> = number of grid points in the 1st coordinate direction\n\ 20 -da_grid_y <ny>, where <ny> = number of grid points in the 2nd coordinate direction\n\ 21 -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise\n\ 22 -lb <value>, lower bound on the variables\n\ 23 -ub <value>, upper bound on the variables\n\n"; 24 25 /* 26 User-defined application context - contains data needed by the 27 application-provided call-back routines, FormJacobian() and 28 FormFunction(). 29 */ 30 31 /* 32 This is a new version of the ../tests/ex8.c code 33 34 Run, for example, with the options ./ex58 -snes_vi_monitor -ksp_monitor -mg_levels_ksp_monitor -pc_type mg -pc_mg_levels 2 -pc_mg_galerkin pmat -ksp_type fgmres 35 36 Or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of 37 multigrid levels, it will be determined automatically based on the number of refinements done) 38 39 ./ex58 -pc_type mg -ksp_monitor -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3 40 -mg_levels_ksp_monitor -snes_vi_monitor -mg_levels_pc_type sor -pc_mg_type full 41 42 */ 43 44 typedef struct { 45 PetscScalar *bottom, *top, *left, *right; 46 PetscScalar lb,ub; 47 } AppCtx; 48 49 /* -------- User-defined Routines --------- */ 50 51 extern PetscErrorCode FormBoundaryConditions(SNES,AppCtx**); 52 extern PetscErrorCode DestroyBoundaryConditions(AppCtx**); 53 extern PetscErrorCode ComputeInitialGuess(SNES,Vec,void*); 54 extern PetscErrorCode FormGradient(SNES,Vec,Vec,void*); 55 extern PetscErrorCode FormJacobian(SNES,Vec,Mat,Mat,void*); 56 extern PetscErrorCode FormBounds(SNES,Vec,Vec); 57 58 int main(int argc, char **argv) 59 { 60 Vec x,r; /* solution and residual vectors */ 61 SNES snes; /* nonlinear solver context */ 62 Mat J; /* Jacobian matrix */ 63 DM da; 64 65 CHKERRQ(PetscInitialize(&argc, &argv, (char*)0, help)); 66 67 /* Create distributed array to manage the 2d grid */ 68 CHKERRQ(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE,DMDA_STENCIL_BOX,4,4,PETSC_DECIDE,PETSC_DECIDE,1,1,NULL,NULL,&da)); 69 CHKERRQ(DMSetFromOptions(da)); 70 CHKERRQ(DMSetUp(da)); 71 72 /* Extract global vectors from DMDA; */ 73 CHKERRQ(DMCreateGlobalVector(da,&x)); 74 CHKERRQ(VecDuplicate(x, &r)); 75 76 CHKERRQ(DMSetMatType(da,MATAIJ)); 77 CHKERRQ(DMCreateMatrix(da,&J)); 78 79 /* Create nonlinear solver context */ 80 CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes)); 81 CHKERRQ(SNESSetDM(snes,da)); 82 83 /* Set function evaluation and Jacobian evaluation routines */ 84 CHKERRQ(SNESSetFunction(snes,r,FormGradient,NULL)); 85 CHKERRQ(SNESSetJacobian(snes,J,J,FormJacobian,NULL)); 86 87 CHKERRQ(SNESSetComputeApplicationContext(snes,(PetscErrorCode (*)(SNES,void**))FormBoundaryConditions,(PetscErrorCode (*)(void**))DestroyBoundaryConditions)); 88 89 CHKERRQ(SNESSetComputeInitialGuess(snes,ComputeInitialGuess,NULL)); 90 91 CHKERRQ(SNESVISetComputeVariableBounds(snes,FormBounds)); 92 93 CHKERRQ(SNESSetFromOptions(snes)); 94 95 /* Solve the application */ 96 CHKERRQ(SNESSolve(snes,NULL,x)); 97 98 /* Free memory */ 99 CHKERRQ(VecDestroy(&x)); 100 CHKERRQ(VecDestroy(&r)); 101 CHKERRQ(MatDestroy(&J)); 102 CHKERRQ(SNESDestroy(&snes)); 103 104 /* Free user-created data structures */ 105 CHKERRQ(DMDestroy(&da)); 106 107 CHKERRQ(PetscFinalize()); 108 return 0; 109 } 110 111 /* -------------------------------------------------------------------- */ 112 113 /* FormBounds - sets the upper and lower bounds 114 115 Input Parameters: 116 . snes - the SNES context 117 118 Output Parameters: 119 . xl - lower bounds 120 . xu - upper bounds 121 */ 122 PetscErrorCode FormBounds(SNES snes, Vec xl, Vec xu) 123 { 124 AppCtx *ctx; 125 126 PetscFunctionBeginUser; 127 CHKERRQ(SNESGetApplicationContext(snes,&ctx)); 128 CHKERRQ(VecSet(xl,ctx->lb)); 129 CHKERRQ(VecSet(xu,ctx->ub)); 130 PetscFunctionReturn(0); 131 } 132 133 /* -------------------------------------------------------------------- */ 134 135 /* FormGradient - Evaluates gradient of f. 136 137 Input Parameters: 138 . snes - the SNES context 139 . X - input vector 140 . ptr - optional user-defined context, as set by SNESSetFunction() 141 142 Output Parameters: 143 . G - vector containing the newly evaluated gradient 144 */ 145 PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr) 146 { 147 AppCtx *user; 148 PetscInt i,j; 149 PetscInt mx, my; 150 PetscScalar hx,hy, hydhx, hxdhy; 151 PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 152 PetscScalar df1dxc,df2dxc,df3dxc,df4dxc,df5dxc,df6dxc; 153 PetscScalar **g, **x; 154 PetscInt xs,xm,ys,ym; 155 Vec localX; 156 DM da; 157 158 PetscFunctionBeginUser; 159 CHKERRQ(SNESGetDM(snes,&da)); 160 CHKERRQ(SNESGetApplicationContext(snes,&user)); 161 CHKERRQ(DMDAGetInfo(da,PETSC_IGNORE,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE)); 162 hx = 1.0/(mx+1);hy=1.0/(my+1); hydhx=hy/hx; hxdhy=hx/hy; 163 164 CHKERRQ(VecSet(G,0.0)); 165 166 /* Get local vector */ 167 CHKERRQ(DMGetLocalVector(da,&localX)); 168 /* Get ghost points */ 169 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 170 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 171 /* Get pointer to local vector data */ 172 CHKERRQ(DMDAVecGetArray(da,localX, &x)); 173 CHKERRQ(DMDAVecGetArray(da,G, &g)); 174 175 CHKERRQ(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 176 /* Compute function over the locally owned part of the mesh */ 177 for (j=ys; j < ys+ym; j++) { 178 for (i=xs; i< xs+xm; i++) { 179 180 xc = x[j][i]; 181 xlt=xrb=xl=xr=xb=xt=xc; 182 183 if (i==0) { /* left side */ 184 xl = user->left[j+1]; 185 xlt = user->left[j+2]; 186 } else xl = x[j][i-1]; 187 188 if (j==0) { /* bottom side */ 189 xb = user->bottom[i+1]; 190 xrb = user->bottom[i+2]; 191 } else xb = x[j-1][i]; 192 193 if (i+1 == mx) { /* right side */ 194 xr = user->right[j+1]; 195 xrb = user->right[j]; 196 } else xr = x[j][i+1]; 197 198 if (j+1==0+my) { /* top side */ 199 xt = user->top[i+1]; 200 xlt = user->top[i]; 201 } else xt = x[j+1][i]; 202 203 if (i>0 && j+1<my) xlt = x[j+1][i-1]; /* left top side */ 204 if (j>0 && i+1<mx) xrb = x[j-1][i+1]; /* right bottom */ 205 206 d1 = (xc-xl); 207 d2 = (xc-xr); 208 d3 = (xc-xt); 209 d4 = (xc-xb); 210 d5 = (xr-xrb); 211 d6 = (xrb-xb); 212 d7 = (xlt-xl); 213 d8 = (xt-xlt); 214 215 df1dxc = d1*hydhx; 216 df2dxc = (d1*hydhx + d4*hxdhy); 217 df3dxc = d3*hxdhy; 218 df4dxc = (d2*hydhx + d3*hxdhy); 219 df5dxc = d2*hydhx; 220 df6dxc = d4*hxdhy; 221 222 d1 /= hx; 223 d2 /= hx; 224 d3 /= hy; 225 d4 /= hy; 226 d5 /= hy; 227 d6 /= hx; 228 d7 /= hy; 229 d8 /= hx; 230 231 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7); 232 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4); 233 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8); 234 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2); 235 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5); 236 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6); 237 238 df1dxc /= f1; 239 df2dxc /= f2; 240 df3dxc /= f3; 241 df4dxc /= f4; 242 df5dxc /= f5; 243 df6dxc /= f6; 244 245 g[j][i] = (df1dxc+df2dxc+df3dxc+df4dxc+df5dxc+df6dxc)/2.0; 246 247 } 248 } 249 250 /* Restore vectors */ 251 CHKERRQ(DMDAVecRestoreArray(da,localX, &x)); 252 CHKERRQ(DMDAVecRestoreArray(da,G, &g)); 253 CHKERRQ(DMRestoreLocalVector(da,&localX)); 254 CHKERRQ(PetscLogFlops(67.0*mx*my)); 255 PetscFunctionReturn(0); 256 } 257 258 /* ------------------------------------------------------------------- */ 259 /* 260 FormJacobian - Evaluates Jacobian matrix. 261 262 Input Parameters: 263 . snes - SNES context 264 . X - input vector 265 . ptr - optional user-defined context, as set by SNESSetJacobian() 266 267 Output Parameters: 268 . tH - Jacobian matrix 269 270 */ 271 PetscErrorCode FormJacobian(SNES snes, Vec X, Mat H, Mat tHPre, void *ptr) 272 { 273 AppCtx *user; 274 PetscInt i,j,k; 275 PetscInt mx, my; 276 MatStencil row,col[7]; 277 PetscScalar hx, hy, hydhx, hxdhy; 278 PetscScalar f1,f2,f3,f4,f5,f6,d1,d2,d3,d4,d5,d6,d7,d8,xc,xl,xr,xt,xb,xlt,xrb; 279 PetscScalar hl,hr,ht,hb,hc,htl,hbr; 280 PetscScalar **x, v[7]; 281 PetscBool assembled; 282 PetscInt xs,xm,ys,ym; 283 Vec localX; 284 DM da; 285 286 PetscFunctionBeginUser; 287 CHKERRQ(SNESGetDM(snes,&da)); 288 CHKERRQ(SNESGetApplicationContext(snes,&user)); 289 CHKERRQ(DMDAGetInfo(da,PETSC_IGNORE,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE)); 290 hx = 1.0/(mx+1); hy=1.0/(my+1); hydhx=hy/hx; hxdhy=hx/hy; 291 292 /* Set various matrix options */ 293 CHKERRQ(MatAssembled(H,&assembled)); 294 if (assembled) CHKERRQ(MatZeroEntries(H)); 295 296 /* Get local vector */ 297 CHKERRQ(DMGetLocalVector(da,&localX)); 298 /* Get ghost points */ 299 CHKERRQ(DMGlobalToLocalBegin(da,X,INSERT_VALUES,localX)); 300 CHKERRQ(DMGlobalToLocalEnd(da,X,INSERT_VALUES,localX)); 301 302 /* Get pointers to vector data */ 303 CHKERRQ(DMDAVecGetArray(da,localX, &x)); 304 305 CHKERRQ(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 306 /* Compute Jacobian over the locally owned part of the mesh */ 307 for (j=ys; j< ys+ym; j++) { 308 for (i=xs; i< xs+xm; i++) { 309 xc = x[j][i]; 310 xlt=xrb=xl=xr=xb=xt=xc; 311 312 /* Left */ 313 if (i==0) { 314 xl = user->left[j+1]; 315 xlt = user->left[j+2]; 316 } else xl = x[j][i-1]; 317 318 /* Bottom */ 319 if (j==0) { 320 xb =user->bottom[i+1]; 321 xrb = user->bottom[i+2]; 322 } else xb = x[j-1][i]; 323 324 /* Right */ 325 if (i+1 == mx) { 326 xr =user->right[j+1]; 327 xrb = user->right[j]; 328 } else xr = x[j][i+1]; 329 330 /* Top */ 331 if (j+1==my) { 332 xt =user->top[i+1]; 333 xlt = user->top[i]; 334 } else xt = x[j+1][i]; 335 336 /* Top left */ 337 if (i>0 && j+1<my) xlt = x[j+1][i-1]; 338 339 /* Bottom right */ 340 if (j>0 && i+1<mx) xrb = x[j-1][i+1]; 341 342 d1 = (xc-xl)/hx; 343 d2 = (xc-xr)/hx; 344 d3 = (xc-xt)/hy; 345 d4 = (xc-xb)/hy; 346 d5 = (xrb-xr)/hy; 347 d6 = (xrb-xb)/hx; 348 d7 = (xlt-xl)/hy; 349 d8 = (xlt-xt)/hx; 350 351 f1 = PetscSqrtScalar(1.0 + d1*d1 + d7*d7); 352 f2 = PetscSqrtScalar(1.0 + d1*d1 + d4*d4); 353 f3 = PetscSqrtScalar(1.0 + d3*d3 + d8*d8); 354 f4 = PetscSqrtScalar(1.0 + d3*d3 + d2*d2); 355 f5 = PetscSqrtScalar(1.0 + d2*d2 + d5*d5); 356 f6 = PetscSqrtScalar(1.0 + d4*d4 + d6*d6); 357 358 hl = (-hydhx*(1.0+d7*d7)+d1*d7)/(f1*f1*f1)+ 359 (-hydhx*(1.0+d4*d4)+d1*d4)/(f2*f2*f2); 360 hr = (-hydhx*(1.0+d5*d5)+d2*d5)/(f5*f5*f5)+ 361 (-hydhx*(1.0+d3*d3)+d2*d3)/(f4*f4*f4); 362 ht = (-hxdhy*(1.0+d8*d8)+d3*d8)/(f3*f3*f3)+ 363 (-hxdhy*(1.0+d2*d2)+d2*d3)/(f4*f4*f4); 364 hb = (-hxdhy*(1.0+d6*d6)+d4*d6)/(f6*f6*f6)+ 365 (-hxdhy*(1.0+d1*d1)+d1*d4)/(f2*f2*f2); 366 367 hbr = -d2*d5/(f5*f5*f5) - d4*d6/(f6*f6*f6); 368 htl = -d1*d7/(f1*f1*f1) - d3*d8/(f3*f3*f3); 369 370 hc = hydhx*(1.0+d7*d7)/(f1*f1*f1) + hxdhy*(1.0+d8*d8)/(f3*f3*f3) + 371 hydhx*(1.0+d5*d5)/(f5*f5*f5) + hxdhy*(1.0+d6*d6)/(f6*f6*f6) + 372 (hxdhy*(1.0+d1*d1)+hydhx*(1.0+d4*d4)-2.0*d1*d4)/(f2*f2*f2) + 373 (hxdhy*(1.0+d2*d2)+hydhx*(1.0+d3*d3)-2.0*d2*d3)/(f4*f4*f4); 374 375 hl/=2.0; hr/=2.0; ht/=2.0; hb/=2.0; hbr/=2.0; htl/=2.0; hc/=2.0; 376 377 k =0; 378 row.i = i;row.j= j; 379 /* Bottom */ 380 if (j>0) { 381 v[k] =hb; 382 col[k].i = i; col[k].j=j-1; k++; 383 } 384 385 /* Bottom right */ 386 if (j>0 && i < mx -1) { 387 v[k] =hbr; 388 col[k].i = i+1; col[k].j = j-1; k++; 389 } 390 391 /* left */ 392 if (i>0) { 393 v[k] = hl; 394 col[k].i = i-1; col[k].j = j; k++; 395 } 396 397 /* Centre */ 398 v[k]= hc; col[k].i= row.i; col[k].j = row.j; k++; 399 400 /* Right */ 401 if (i < mx-1) { 402 v[k] = hr; 403 col[k].i= i+1; col[k].j = j;k++; 404 } 405 406 /* Top left */ 407 if (i>0 && j < my-1) { 408 v[k] = htl; 409 col[k].i = i-1;col[k].j = j+1; k++; 410 } 411 412 /* Top */ 413 if (j < my-1) { 414 v[k] = ht; 415 col[k].i = i; col[k].j = j+1; k++; 416 } 417 418 CHKERRQ(MatSetValuesStencil(H,1,&row,k,col,v,INSERT_VALUES)); 419 } 420 } 421 422 /* Assemble the matrix */ 423 CHKERRQ(MatAssemblyBegin(H,MAT_FINAL_ASSEMBLY)); 424 CHKERRQ(DMDAVecRestoreArray(da,localX,&x)); 425 CHKERRQ(MatAssemblyEnd(H,MAT_FINAL_ASSEMBLY)); 426 CHKERRQ(DMRestoreLocalVector(da,&localX)); 427 428 CHKERRQ(PetscLogFlops(199.0*mx*my)); 429 PetscFunctionReturn(0); 430 } 431 432 /* ------------------------------------------------------------------- */ 433 /* 434 FormBoundaryConditions - Calculates the boundary conditions for 435 the region. 436 437 Input Parameter: 438 . user - user-defined application context 439 440 Output Parameter: 441 . user - user-defined application context 442 */ 443 PetscErrorCode FormBoundaryConditions(SNES snes,AppCtx **ouser) 444 { 445 PetscInt i,j,k,limit=0,maxits=5; 446 PetscInt mx,my; 447 PetscInt bsize=0, lsize=0, tsize=0, rsize=0; 448 PetscScalar one =1.0, two=2.0, three=3.0; 449 PetscScalar det,hx,hy,xt=0,yt=0; 450 PetscReal fnorm, tol=1e-10; 451 PetscScalar u1,u2,nf1,nf2,njac11,njac12,njac21,njac22; 452 PetscScalar b=-0.5, t=0.5, l=-0.5, r=0.5; 453 PetscScalar *boundary; 454 AppCtx *user; 455 DM da; 456 457 PetscFunctionBeginUser; 458 CHKERRQ(SNESGetDM(snes,&da)); 459 CHKERRQ(PetscNew(&user)); 460 *ouser = user; 461 user->lb = .05; 462 user->ub = PETSC_INFINITY; 463 CHKERRQ(DMDAGetInfo(da,PETSC_IGNORE,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE)); 464 465 /* Check if lower and upper bounds are set */ 466 CHKERRQ(PetscOptionsGetScalar(NULL,NULL, "-lb", &user->lb, 0)); 467 CHKERRQ(PetscOptionsGetScalar(NULL,NULL, "-ub", &user->ub, 0)); 468 bsize=mx+2; lsize=my+2; rsize=my+2; tsize=mx+2; 469 470 CHKERRQ(PetscMalloc1(bsize, &user->bottom)); 471 CHKERRQ(PetscMalloc1(tsize, &user->top)); 472 CHKERRQ(PetscMalloc1(lsize, &user->left)); 473 CHKERRQ(PetscMalloc1(rsize, &user->right)); 474 475 hx= (r-l)/(mx+1.0); hy=(t-b)/(my+1.0); 476 477 for (j=0; j<4; j++) { 478 if (j==0) { 479 yt = b; 480 xt = l; 481 limit = bsize; 482 boundary = user->bottom; 483 } else if (j==1) { 484 yt = t; 485 xt = l; 486 limit = tsize; 487 boundary = user->top; 488 } else if (j==2) { 489 yt = b; 490 xt = l; 491 limit = lsize; 492 boundary = user->left; 493 } else { /* if (j==3) */ 494 yt = b; 495 xt = r; 496 limit = rsize; 497 boundary = user->right; 498 } 499 500 for (i=0; i<limit; i++) { 501 u1=xt; 502 u2=-yt; 503 for (k=0; k<maxits; k++) { 504 nf1 = u1 + u1*u2*u2 - u1*u1*u1/three-xt; 505 nf2 = -u2 - u1*u1*u2 + u2*u2*u2/three-yt; 506 fnorm = PetscRealPart(PetscSqrtScalar(nf1*nf1+nf2*nf2)); 507 if (fnorm <= tol) break; 508 njac11=one+u2*u2-u1*u1; 509 njac12=two*u1*u2; 510 njac21=-two*u1*u2; 511 njac22=-one - u1*u1 + u2*u2; 512 det = njac11*njac22-njac21*njac12; 513 u1 = u1-(njac22*nf1-njac12*nf2)/det; 514 u2 = u2-(njac11*nf2-njac21*nf1)/det; 515 } 516 517 boundary[i]=u1*u1-u2*u2; 518 if (j==0 || j==1) xt=xt+hx; 519 else yt=yt+hy; /* if (j==2 || j==3) */ 520 } 521 } 522 PetscFunctionReturn(0); 523 } 524 525 PetscErrorCode DestroyBoundaryConditions(AppCtx **ouser) 526 { 527 AppCtx *user = *ouser; 528 529 PetscFunctionBeginUser; 530 CHKERRQ(PetscFree(user->bottom)); 531 CHKERRQ(PetscFree(user->top)); 532 CHKERRQ(PetscFree(user->left)); 533 CHKERRQ(PetscFree(user->right)); 534 CHKERRQ(PetscFree(*ouser)); 535 PetscFunctionReturn(0); 536 } 537 538 /* ------------------------------------------------------------------- */ 539 /* 540 ComputeInitialGuess - Calculates the initial guess 541 542 Input Parameters: 543 . user - user-defined application context 544 . X - vector for initial guess 545 546 Output Parameters: 547 . X - newly computed initial guess 548 */ 549 PetscErrorCode ComputeInitialGuess(SNES snes, Vec X,void *dummy) 550 { 551 PetscInt i,j,mx,my; 552 DM da; 553 AppCtx *user; 554 PetscScalar **x; 555 PetscInt xs,xm,ys,ym; 556 557 PetscFunctionBeginUser; 558 CHKERRQ(SNESGetDM(snes,&da)); 559 CHKERRQ(SNESGetApplicationContext(snes,&user)); 560 561 CHKERRQ(DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL)); 562 CHKERRQ(DMDAGetInfo(da,PETSC_IGNORE,&mx,&my,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE)); 563 564 /* Get pointers to vector data */ 565 CHKERRQ(DMDAVecGetArray(da,X,&x)); 566 /* Perform local computations */ 567 for (j=ys; j<ys+ym; j++) { 568 for (i=xs; i< xs+xm; i++) { 569 x[j][i] = (((j+1.0)*user->bottom[i+1]+(my-j+1.0)*user->top[i+1])/(my+2.0)+((i+1.0)*user->left[j+1]+(mx-i+1.0)*user->right[j+1])/(mx+2.0))/2.0; 570 } 571 } 572 /* Restore vectors */ 573 CHKERRQ(DMDAVecRestoreArray(da,X,&x)); 574 PetscFunctionReturn(0); 575 } 576 577 /*TEST 578 579 test: 580 args: -snes_type vinewtonrsls -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason 581 requires: !single 582 583 test: 584 suffix: 2 585 args: -snes_type vinewtonssls -pc_type mg -ksp_monitor_short -pc_mg_galerkin pmat -da_refine 5 -snes_vi_monitor -pc_mg_type full -snes_max_it 100 -snes_converged_reason 586 requires: !single 587 588 TEST*/ 589