Actual source code: ex58.c


  2: #include <petscsnes.h>
  3: #include <petscdm.h>
  4: #include <petscdmda.h>

  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";

 25: /*
 26:    User-defined application context - contains data needed by the
 27:    application-provided call-back routines, FormJacobian() and
 28:    FormFunction().
 29: */

 31: /*
 32:      This is a new version of the ../tests/ex8.c code

 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

 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)

 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

 42: */

 44: typedef struct {
 45:   PetscScalar *bottom, *top, *left, *right;
 46:   PetscScalar  lb, ub;
 47: } AppCtx;

 49: /* -------- User-defined Routines --------- */

 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);

 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;

 66:   PetscInitialize(&argc, &argv, (char *)0, help);

 68:   /* Create distributed array to manage the 2d grid */
 69:   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:   DMSetFromOptions(da);
 71:   DMSetUp(da);

 73:   /* Extract global vectors from DMDA; */
 74:   DMCreateGlobalVector(da, &x);
 75:   VecDuplicate(x, &r);

 77:   DMSetMatType(da, MATAIJ);
 78:   DMCreateMatrix(da, &J);

 80:   /* Create nonlinear solver context */
 81:   SNESCreate(PETSC_COMM_WORLD, &snes);
 82:   SNESSetDM(snes, da);

 84:   /*  Set function evaluation and Jacobian evaluation  routines */
 85:   SNESSetFunction(snes, r, FormGradient, NULL);
 86:   SNESSetJacobian(snes, J, J, FormJacobian, NULL);

 88:   SNESSetComputeApplicationContext(snes, (PetscErrorCode(*)(SNES, void **))FormBoundaryConditions, (PetscErrorCode(*)(void **))DestroyBoundaryConditions);

 90:   SNESSetComputeInitialGuess(snes, ComputeInitialGuess, NULL);

 92:   SNESVISetComputeVariableBounds(snes, FormBounds);

 94:   SNESSetFromOptions(snes);

 96:   /* Solve the application */
 97:   SNESSolve(snes, NULL, x);

 99:   /* Free memory */
100:   VecDestroy(&x);
101:   VecDestroy(&r);
102:   MatDestroy(&J);
103:   SNESDestroy(&snes);

105:   /* Free user-created data structures */
106:   DMDestroy(&da);

108:   PetscFinalize();
109:   return 0;
110: }

112: /* -------------------------------------------------------------------- */

114: /*  FormBounds - sets the upper and lower bounds

116:     Input Parameters:
117: .   snes  - the SNES context

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;

128:   SNESGetApplicationContext(snes, &ctx);
129:   VecSet(xl, ctx->lb);
130:   VecSet(xu, ctx->ub);
131:   return 0;
132: }

134: /* -------------------------------------------------------------------- */

136: /*  FormGradient - Evaluates gradient of f.

138:     Input Parameters:
139: .   snes  - the SNES context
140: .   X     - input vector
141: .   ptr   - optional user-defined context, as set by SNESSetFunction()

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;

160:   SNESGetDM(snes, &da);
161:   SNESGetApplicationContext(snes, &user);
162:   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);
164:   hy    = 1.0 / (my + 1);
165:   hydhx = hy / hx;
166:   hxdhy = hx / hy;

168:   VecSet(G, 0.0);

170:   /* Get local vector */
171:   DMGetLocalVector(da, &localX);
172:   /* Get ghost points */
173:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
174:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);
175:   /* Get pointer to local vector data */
176:   DMDAVecGetArray(da, localX, &x);
177:   DMDAVecGetArray(da, G, &g);

179:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
180:   /* Compute function over the locally owned part of the mesh */
181:   for (j = ys; j < ys + ym; j++) {
182:     for (i = xs; i < xs + xm; i++) {
183:       xc  = x[j][i];
184:       xlt = xrb = xl = xr = xb = xt = xc;

186:       if (i == 0) { /* left side */
187:         xl  = user->left[j + 1];
188:         xlt = user->left[j + 2];
189:       } else xl = x[j][i - 1];

191:       if (j == 0) { /* bottom side */
192:         xb  = user->bottom[i + 1];
193:         xrb = user->bottom[i + 2];
194:       } else xb = x[j - 1][i];

196:       if (i + 1 == mx) { /* right side */
197:         xr  = user->right[j + 1];
198:         xrb = user->right[j];
199:       } else xr = x[j][i + 1];

201:       if (j + 1 == 0 + my) { /* top side */
202:         xt  = user->top[i + 1];
203:         xlt = user->top[i];
204:       } else xt = x[j + 1][i];

206:       if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1]; /* left top side */
207:       if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1]; /* right bottom */

209:       d1 = (xc - xl);
210:       d2 = (xc - xr);
211:       d3 = (xc - xt);
212:       d4 = (xc - xb);
213:       d5 = (xr - xrb);
214:       d6 = (xrb - xb);
215:       d7 = (xlt - xl);
216:       d8 = (xt - xlt);

218:       df1dxc = d1 * hydhx;
219:       df2dxc = (d1 * hydhx + d4 * hxdhy);
220:       df3dxc = d3 * hxdhy;
221:       df4dxc = (d2 * hydhx + d3 * hxdhy);
222:       df5dxc = d2 * hydhx;
223:       df6dxc = d4 * hxdhy;

225:       d1 /= hx;
226:       d2 /= hx;
227:       d3 /= hy;
228:       d4 /= hy;
229:       d5 /= hy;
230:       d6 /= hx;
231:       d7 /= hy;
232:       d8 /= hx;

234:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
235:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
236:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
237:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
238:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
239:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

241:       df1dxc /= f1;
242:       df2dxc /= f2;
243:       df3dxc /= f3;
244:       df4dxc /= f4;
245:       df5dxc /= f5;
246:       df6dxc /= f6;

248:       g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
249:     }
250:   }

252:   /* Restore vectors */
253:   DMDAVecRestoreArray(da, localX, &x);
254:   DMDAVecRestoreArray(da, G, &g);
255:   DMRestoreLocalVector(da, &localX);
256:   PetscLogFlops(67.0 * mx * my);
257:   return 0;
258: }

260: /* ------------------------------------------------------------------- */
261: /*
262:    FormJacobian - Evaluates Jacobian matrix.

264:    Input Parameters:
265: .  snes - SNES context
266: .  X    - input vector
267: .  ptr  - optional user-defined context, as set by SNESSetJacobian()

269:    Output Parameters:
270: .  tH    - Jacobian matrix

272: */
273: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat H, Mat tHPre, void *ptr)
274: {
275:   AppCtx       *user;
276:   PetscInt      i, j, k;
277:   PetscInt      mx, my;
278:   MatStencil    row, col[7];
279:   PetscScalar   hx, hy, hydhx, hxdhy;
280:   PetscScalar   f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
281:   PetscScalar   hl, hr, ht, hb, hc, htl, hbr;
282:   PetscScalar **x, v[7];
283:   PetscBool     assembled;
284:   PetscInt      xs, xm, ys, ym;
285:   Vec           localX;
286:   DM            da;

289:   SNESGetDM(snes, &da);
290:   SNESGetApplicationContext(snes, &user);
291:   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);
292:   hx    = 1.0 / (mx + 1);
293:   hy    = 1.0 / (my + 1);
294:   hydhx = hy / hx;
295:   hxdhy = hx / hy;

297:   /* Set various matrix options */
298:   MatAssembled(H, &assembled);
299:   if (assembled) MatZeroEntries(H);

301:   /* Get local vector */
302:   DMGetLocalVector(da, &localX);
303:   /* Get ghost points */
304:   DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX);
305:   DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX);

307:   /* Get pointers to vector data */
308:   DMDAVecGetArray(da, localX, &x);

310:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
311:   /* Compute Jacobian over the locally owned part of the mesh */
312:   for (j = ys; j < ys + ym; j++) {
313:     for (i = xs; i < xs + xm; i++) {
314:       xc  = x[j][i];
315:       xlt = xrb = xl = xr = xb = xt = xc;

317:       /* Left */
318:       if (i == 0) {
319:         xl  = user->left[j + 1];
320:         xlt = user->left[j + 2];
321:       } else xl = x[j][i - 1];

323:       /* Bottom */
324:       if (j == 0) {
325:         xb  = user->bottom[i + 1];
326:         xrb = user->bottom[i + 2];
327:       } else xb = x[j - 1][i];

329:       /* Right */
330:       if (i + 1 == mx) {
331:         xr  = user->right[j + 1];
332:         xrb = user->right[j];
333:       } else xr = x[j][i + 1];

335:       /* Top */
336:       if (j + 1 == my) {
337:         xt  = user->top[i + 1];
338:         xlt = user->top[i];
339:       } else xt = x[j + 1][i];

341:       /* Top left */
342:       if (i > 0 && j + 1 < my) xlt = x[j + 1][i - 1];

344:       /* Bottom right */
345:       if (j > 0 && i + 1 < mx) xrb = x[j - 1][i + 1];

347:       d1 = (xc - xl) / hx;
348:       d2 = (xc - xr) / hx;
349:       d3 = (xc - xt) / hy;
350:       d4 = (xc - xb) / hy;
351:       d5 = (xrb - xr) / hy;
352:       d6 = (xrb - xb) / hx;
353:       d7 = (xlt - xl) / hy;
354:       d8 = (xlt - xt) / hx;

356:       f1 = PetscSqrtScalar(1.0 + d1 * d1 + d7 * d7);
357:       f2 = PetscSqrtScalar(1.0 + d1 * d1 + d4 * d4);
358:       f3 = PetscSqrtScalar(1.0 + d3 * d3 + d8 * d8);
359:       f4 = PetscSqrtScalar(1.0 + d3 * d3 + d2 * d2);
360:       f5 = PetscSqrtScalar(1.0 + d2 * d2 + d5 * d5);
361:       f6 = PetscSqrtScalar(1.0 + d4 * d4 + d6 * d6);

363:       hl = (-hydhx * (1.0 + d7 * d7) + d1 * d7) / (f1 * f1 * f1) + (-hydhx * (1.0 + d4 * d4) + d1 * d4) / (f2 * f2 * f2);
364:       hr = (-hydhx * (1.0 + d5 * d5) + d2 * d5) / (f5 * f5 * f5) + (-hydhx * (1.0 + d3 * d3) + d2 * d3) / (f4 * f4 * f4);
365:       ht = (-hxdhy * (1.0 + d8 * d8) + d3 * d8) / (f3 * f3 * f3) + (-hxdhy * (1.0 + d2 * d2) + d2 * d3) / (f4 * f4 * f4);
366:       hb = (-hxdhy * (1.0 + d6 * d6) + d4 * d6) / (f6 * f6 * f6) + (-hxdhy * (1.0 + d1 * d1) + d1 * d4) / (f2 * f2 * f2);

368:       hbr = -d2 * d5 / (f5 * f5 * f5) - d4 * d6 / (f6 * f6 * f6);
369:       htl = -d1 * d7 / (f1 * f1 * f1) - d3 * d8 / (f3 * f3 * f3);

371:       hc = hydhx * (1.0 + d7 * d7) / (f1 * f1 * f1) + hxdhy * (1.0 + d8 * d8) / (f3 * f3 * f3) + hydhx * (1.0 + d5 * d5) / (f5 * f5 * f5) + hxdhy * (1.0 + d6 * d6) / (f6 * f6 * f6) + (hxdhy * (1.0 + d1 * d1) + hydhx * (1.0 + d4 * d4) - 2.0 * d1 * d4) / (f2 * f2 * f2) + (hxdhy * (1.0 + d2 * d2) + hydhx * (1.0 + d3 * d3) - 2.0 * d2 * d3) / (f4 * f4 * f4);

373:       hl /= 2.0;
374:       hr /= 2.0;
375:       ht /= 2.0;
376:       hb /= 2.0;
377:       hbr /= 2.0;
378:       htl /= 2.0;
379:       hc /= 2.0;

381:       k     = 0;
382:       row.i = i;
383:       row.j = j;
384:       /* Bottom */
385:       if (j > 0) {
386:         v[k]     = hb;
387:         col[k].i = i;
388:         col[k].j = j - 1;
389:         k++;
390:       }

392:       /* Bottom right */
393:       if (j > 0 && i < mx - 1) {
394:         v[k]     = hbr;
395:         col[k].i = i + 1;
396:         col[k].j = j - 1;
397:         k++;
398:       }

400:       /* left */
401:       if (i > 0) {
402:         v[k]     = hl;
403:         col[k].i = i - 1;
404:         col[k].j = j;
405:         k++;
406:       }

408:       /* Centre */
409:       v[k]     = hc;
410:       col[k].i = row.i;
411:       col[k].j = row.j;
412:       k++;

414:       /* Right */
415:       if (i < mx - 1) {
416:         v[k]     = hr;
417:         col[k].i = i + 1;
418:         col[k].j = j;
419:         k++;
420:       }

422:       /* Top left */
423:       if (i > 0 && j < my - 1) {
424:         v[k]     = htl;
425:         col[k].i = i - 1;
426:         col[k].j = j + 1;
427:         k++;
428:       }

430:       /* Top */
431:       if (j < my - 1) {
432:         v[k]     = ht;
433:         col[k].i = i;
434:         col[k].j = j + 1;
435:         k++;
436:       }

438:       MatSetValuesStencil(H, 1, &row, k, col, v, INSERT_VALUES);
439:     }
440:   }

442:   /* Assemble the matrix */
443:   MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY);
444:   DMDAVecRestoreArray(da, localX, &x);
445:   MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY);
446:   DMRestoreLocalVector(da, &localX);

448:   PetscLogFlops(199.0 * mx * my);
449:   return 0;
450: }

452: /* ------------------------------------------------------------------- */
453: /*
454:    FormBoundaryConditions -  Calculates the boundary conditions for
455:    the region.

457:    Input Parameter:
458: .  user - user-defined application context

460:    Output Parameter:
461: .  user - user-defined application context
462: */
463: PetscErrorCode FormBoundaryConditions(SNES snes, AppCtx **ouser)
464: {
465:   PetscInt     i, j, k, limit = 0, maxits = 5;
466:   PetscInt     mx, my;
467:   PetscInt     bsize = 0, lsize = 0, tsize = 0, rsize = 0;
468:   PetscScalar  one = 1.0, two = 2.0, three = 3.0;
469:   PetscScalar  det, hx, hy, xt = 0, yt = 0;
470:   PetscReal    fnorm, tol = 1e-10;
471:   PetscScalar  u1, u2, nf1, nf2, njac11, njac12, njac21, njac22;
472:   PetscScalar  b = -0.5, t = 0.5, l = -0.5, r = 0.5;
473:   PetscScalar *boundary;
474:   AppCtx      *user;
475:   DM           da;

478:   SNESGetDM(snes, &da);
479:   PetscNew(&user);
480:   *ouser   = user;
481:   user->lb = .05;
482:   user->ub = PETSC_INFINITY;
483:   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);

485:   /* Check if lower and upper bounds are set */
486:   PetscOptionsGetScalar(NULL, NULL, "-lb", &user->lb, 0);
487:   PetscOptionsGetScalar(NULL, NULL, "-ub", &user->ub, 0);
488:   bsize = mx + 2;
489:   lsize = my + 2;
490:   rsize = my + 2;
491:   tsize = mx + 2;

493:   PetscMalloc1(bsize, &user->bottom);
494:   PetscMalloc1(tsize, &user->top);
495:   PetscMalloc1(lsize, &user->left);
496:   PetscMalloc1(rsize, &user->right);

498:   hx = (r - l) / (mx + 1.0);
499:   hy = (t - b) / (my + 1.0);

501:   for (j = 0; j < 4; j++) {
502:     if (j == 0) {
503:       yt       = b;
504:       xt       = l;
505:       limit    = bsize;
506:       boundary = user->bottom;
507:     } else if (j == 1) {
508:       yt       = t;
509:       xt       = l;
510:       limit    = tsize;
511:       boundary = user->top;
512:     } else if (j == 2) {
513:       yt       = b;
514:       xt       = l;
515:       limit    = lsize;
516:       boundary = user->left;
517:     } else { /* if  (j==3) */
518:       yt       = b;
519:       xt       = r;
520:       limit    = rsize;
521:       boundary = user->right;
522:     }

524:     for (i = 0; i < limit; i++) {
525:       u1 = xt;
526:       u2 = -yt;
527:       for (k = 0; k < maxits; k++) {
528:         nf1   = u1 + u1 * u2 * u2 - u1 * u1 * u1 / three - xt;
529:         nf2   = -u2 - u1 * u1 * u2 + u2 * u2 * u2 / three - yt;
530:         fnorm = PetscRealPart(PetscSqrtScalar(nf1 * nf1 + nf2 * nf2));
531:         if (fnorm <= tol) break;
532:         njac11 = one + u2 * u2 - u1 * u1;
533:         njac12 = two * u1 * u2;
534:         njac21 = -two * u1 * u2;
535:         njac22 = -one - u1 * u1 + u2 * u2;
536:         det    = njac11 * njac22 - njac21 * njac12;
537:         u1     = u1 - (njac22 * nf1 - njac12 * nf2) / det;
538:         u2     = u2 - (njac11 * nf2 - njac21 * nf1) / det;
539:       }

541:       boundary[i] = u1 * u1 - u2 * u2;
542:       if (j == 0 || j == 1) xt = xt + hx;
543:       else yt = yt + hy; /* if (j==2 || j==3) */
544:     }
545:   }
546:   return 0;
547: }

549: PetscErrorCode DestroyBoundaryConditions(AppCtx **ouser)
550: {
551:   AppCtx *user = *ouser;

554:   PetscFree(user->bottom);
555:   PetscFree(user->top);
556:   PetscFree(user->left);
557:   PetscFree(user->right);
558:   PetscFree(*ouser);
559:   return 0;
560: }

562: /* ------------------------------------------------------------------- */
563: /*
564:    ComputeInitialGuess - Calculates the initial guess

566:    Input Parameters:
567: .  user - user-defined application context
568: .  X - vector for initial guess

570:    Output Parameters:
571: .  X - newly computed initial guess
572: */
573: PetscErrorCode ComputeInitialGuess(SNES snes, Vec X, void *dummy)
574: {
575:   PetscInt      i, j, mx, my;
576:   DM            da;
577:   AppCtx       *user;
578:   PetscScalar **x;
579:   PetscInt      xs, xm, ys, ym;

582:   SNESGetDM(snes, &da);
583:   SNESGetApplicationContext(snes, &user);

585:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);
586:   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);

588:   /* Get pointers to vector data */
589:   DMDAVecGetArray(da, X, &x);
590:   /* Perform local computations */
591:   for (j = ys; j < ys + ym; j++) {
592:     for (i = xs; i < xs + xm; i++) 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;
593:   }
594:   /* Restore vectors */
595:   DMDAVecRestoreArray(da, X, &x);
596:   return 0;
597: }

599: /*TEST

601:    test:
602:       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
603:       requires: !single

605:    test:
606:       suffix: 2
607:       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
608:       requires: !single

610: TEST*/