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