Actual source code: ex58.c

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

  5: static const char help[] = "Parallel version of the minimum surface area problem in 2D using DMDA.\n\
  6:  It solves a system of nonlinear equations in mixed\n\
  7: complementarity form.This example is based on a\n\
  8: problem from the MINPACK-2 test suite.  Given a rectangular 2-D domain and\n\
  9: boundary values along the edges of the domain, the objective is to find the\n\
 10: surface with the minimal area that satisfies the boundary conditions.\n\
 11: This application solves this problem using complimentarity -- We are actually\n\
 12: solving the system  (grad f)_i >= 0, if x_i == l_i \n\
 13:                     (grad f)_i = 0, if l_i < x_i < u_i \n\
 14:                     (grad f)_i <= 0, if x_i == u_i  \n\
 15: where f is the function to be minimized. \n\
 16: \n\
 17: The command line options are:\n\
 18:   -da_grid_x <nx>, where <nx> = number of grid points in the 1st coordinate direction\n\
 19:   -da_grid_y <ny>, where <ny> = number of grid points in the 2nd coordinate direction\n\
 20:   -start <st>, where <st> =0 for zero vector, and an average of the boundary conditions otherwise\n\
 21:   -lb <value>, lower bound on the variables\n\
 22:   -ub <value>, upper bound on the variables\n\n";

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

 30: /*
 31:      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

 33:      Or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
 34:          multigrid levels, it will be determined automatically based on the number of refinements done)

 36:       ./ex58 -pc_type mg -ksp_monitor  -snes_view -pc_mg_galerkin pmat -snes_grid_sequence 3
 37:              -mg_levels_ksp_monitor -snes_vi_monitor -mg_levels_pc_type sor -pc_mg_type full

 39: */

 41: typedef struct {
 42:   PetscScalar *bottom, *top, *left, *right;
 43:   PetscScalar  lb, ub;
 44: } AppCtx;

 46: /* -------- User-defined Routines --------- */

 48: extern PetscErrorCode FormBoundaryConditions(SNES, AppCtx **);
 49: extern PetscErrorCode DestroyBoundaryConditions(AppCtx **);
 50: extern PetscErrorCode ComputeInitialGuess(SNES, Vec, void *);
 51: extern PetscErrorCode FormGradient(SNES, Vec, Vec, void *);
 52: extern PetscErrorCode FormJacobian(SNES, Vec, Mat, Mat, void *);
 53: extern PetscErrorCode FormBounds(SNES, Vec, Vec);

 55: int main(int argc, char **argv)
 56: {
 57:   Vec  x, r; /* solution and residual vectors */
 58:   SNES snes; /* nonlinear solver context */
 59:   Mat  J;    /* Jacobian matrix */
 60:   DM   da;

 62:   PetscFunctionBeginUser;
 63:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

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

 70:   /* Extract global vectors from DMDA; */
 71:   PetscCall(DMCreateGlobalVector(da, &x));
 72:   PetscCall(VecDuplicate(x, &r));

 74:   PetscCall(DMSetMatType(da, MATAIJ));
 75:   PetscCall(DMCreateMatrix(da, &J));

 77:   /* Create nonlinear solver context */
 78:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
 79:   PetscCall(SNESSetDM(snes, da));

 81:   /*  Set function evaluation and Jacobian evaluation  routines */
 82:   PetscCall(SNESSetFunction(snes, r, FormGradient, NULL));
 83:   PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, NULL));

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

 87:   PetscCall(SNESSetComputeInitialGuess(snes, ComputeInitialGuess, NULL));

 89:   PetscCall(SNESVISetComputeVariableBounds(snes, FormBounds));

 91:   PetscCall(SNESSetFromOptions(snes));

 93:   /* Solve the application */
 94:   PetscCall(SNESSolve(snes, NULL, x));

 96:   /* Free memory */
 97:   PetscCall(VecDestroy(&x));
 98:   PetscCall(VecDestroy(&r));
 99:   PetscCall(MatDestroy(&J));
100:   PetscCall(SNESDestroy(&snes));

102:   /* Free user-created data structures */
103:   PetscCall(DMDestroy(&da));

105:   PetscCall(PetscFinalize());
106:   return 0;
107: }

109: /* -------------------------------------------------------------------- */

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

113:     Input Parameters:
114: .   snes  - the SNES context

116:     Output Parameters:
117: .   xl - lower bounds
118: .   xu - upper bounds
119: */
120: PetscErrorCode FormBounds(SNES snes, Vec xl, Vec xu)
121: {
122:   AppCtx *ctx;

124:   PetscFunctionBeginUser;
125:   PetscCall(SNESGetApplicationContext(snes, &ctx));
126:   PetscCall(VecSet(xl, ctx->lb));
127:   PetscCall(VecSet(xu, ctx->ub));
128:   PetscFunctionReturn(PETSC_SUCCESS);
129: }

131: /* -------------------------------------------------------------------- */

133: /*  FormGradient - Evaluates gradient of f.

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

140:     Output Parameters:
141: .   G - vector containing the newly evaluated gradient
142: */
143: PetscErrorCode FormGradient(SNES snes, Vec X, Vec G, void *ptr)
144: {
145:   AppCtx       *user;
146:   PetscInt      i, j;
147:   PetscInt      mx, my;
148:   PetscScalar   hx, hy, hydhx, hxdhy;
149:   PetscScalar   f1, f2, f3, f4, f5, f6, d1, d2, d3, d4, d5, d6, d7, d8, xc, xl, xr, xt, xb, xlt, xrb;
150:   PetscScalar   df1dxc, df2dxc, df3dxc, df4dxc, df5dxc, df6dxc;
151:   PetscScalar **g, **x;
152:   PetscInt      xs, xm, ys, ym;
153:   Vec           localX;
154:   DM            da;

156:   PetscFunctionBeginUser;
157:   PetscCall(SNESGetDM(snes, &da));
158:   PetscCall(SNESGetApplicationContext(snes, &user));
159:   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));
160:   hx    = 1.0 / (mx + 1);
161:   hy    = 1.0 / (my + 1);
162:   hydhx = hy / hx;
163:   hxdhy = hx / hy;

165:   PetscCall(VecSet(G, 0.0));

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

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:       xc  = x[j][i];
181:       xlt = xrb = xl = xr = xb = xt = xc;

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

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

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

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

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 */

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

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;

222:       d1 /= hx;
223:       d2 /= hx;
224:       d3 /= hy;
225:       d4 /= hy;
226:       d5 /= hy;
227:       d6 /= hx;
228:       d7 /= hy;
229:       d8 /= hx;

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

238:       df1dxc /= f1;
239:       df2dxc /= f2;
240:       df3dxc /= f3;
241:       df4dxc /= f4;
242:       df5dxc /= f5;
243:       df6dxc /= f6;

245:       g[j][i] = (df1dxc + df2dxc + df3dxc + df4dxc + df5dxc + df6dxc) / 2.0;
246:     }
247:   }

249:   /* Restore vectors */
250:   PetscCall(DMDAVecRestoreArray(da, localX, &x));
251:   PetscCall(DMDAVecRestoreArray(da, G, &g));
252:   PetscCall(DMRestoreLocalVector(da, &localX));
253:   PetscCall(PetscLogFlops(67.0 * mx * my));
254:   PetscFunctionReturn(PETSC_SUCCESS);
255: }

257: /* ------------------------------------------------------------------- */
258: /*
259:    FormJacobian - Evaluates Jacobian matrix.

261:    Input Parameters:
262: .  snes - SNES context
263: .  X    - input vector
264: .  ptr  - optional user-defined context, as set by SNESSetJacobian()

266:    Output Parameters:
267: .  tH    - Jacobian matrix

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

285:   PetscFunctionBeginUser;
286:   PetscCall(SNESGetDM(snes, &da));
287:   PetscCall(SNESGetApplicationContext(snes, &user));
288:   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));
289:   hx    = 1.0 / (mx + 1);
290:   hy    = 1.0 / (my + 1);
291:   hydhx = hy / hx;
292:   hxdhy = hx / hy;

294:   /* Set various matrix options */
295:   PetscCall(MatAssembled(H, &assembled));
296:   if (assembled) PetscCall(MatZeroEntries(H));

298:   /* Get local vector */
299:   PetscCall(DMGetLocalVector(da, &localX));
300:   /* Get ghost points */
301:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
302:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

304:   /* Get pointers to vector data */
305:   PetscCall(DMDAVecGetArray(da, localX, &x));

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

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

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

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

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

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

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

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

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

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

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

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

370:       hl /= 2.0;
371:       hr /= 2.0;
372:       ht /= 2.0;
373:       hb /= 2.0;
374:       hbr /= 2.0;
375:       htl /= 2.0;
376:       hc /= 2.0;

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

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

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

405:       /* Centre */
406:       v[k]     = hc;
407:       col[k].i = row.i;
408:       col[k].j = row.j;
409:       k++;

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

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

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

435:       PetscCall(MatSetValuesStencil(H, 1, &row, k, col, v, INSERT_VALUES));
436:     }
437:   }

439:   /* Assemble the matrix */
440:   PetscCall(MatAssemblyBegin(H, MAT_FINAL_ASSEMBLY));
441:   PetscCall(DMDAVecRestoreArray(da, localX, &x));
442:   PetscCall(MatAssemblyEnd(H, MAT_FINAL_ASSEMBLY));
443:   PetscCall(DMRestoreLocalVector(da, &localX));

445:   PetscCall(PetscLogFlops(199.0 * mx * my));
446:   PetscFunctionReturn(PETSC_SUCCESS);
447: }

449: /* ------------------------------------------------------------------- */
450: /*
451:    FormBoundaryConditions -  Calculates the boundary conditions for
452:    the region.

454:    Input Parameter:
455: .  user - user-defined application context

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

474:   PetscFunctionBeginUser;
475:   PetscCall(SNESGetDM(snes, &da));
476:   PetscCall(PetscNew(&user));
477:   *ouser   = user;
478:   user->lb = .05;
479:   user->ub = PETSC_INFINITY;
480:   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));

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

490:   PetscCall(PetscMalloc1(bsize, &user->bottom));
491:   PetscCall(PetscMalloc1(tsize, &user->top));
492:   PetscCall(PetscMalloc1(lsize, &user->left));
493:   PetscCall(PetscMalloc1(rsize, &user->right));

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

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

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

538:       boundary[i] = u1 * u1 - u2 * u2;
539:       if (j == 0 || j == 1) xt = xt + hx;
540:       else yt = yt + hy; /* if (j==2 || j==3) */
541:     }
542:   }
543:   PetscFunctionReturn(PETSC_SUCCESS);
544: }

546: PetscErrorCode DestroyBoundaryConditions(AppCtx **ouser)
547: {
548:   AppCtx *user = *ouser;

550:   PetscFunctionBeginUser;
551:   PetscCall(PetscFree(user->bottom));
552:   PetscCall(PetscFree(user->top));
553:   PetscCall(PetscFree(user->left));
554:   PetscCall(PetscFree(user->right));
555:   PetscCall(PetscFree(*ouser));
556:   PetscFunctionReturn(PETSC_SUCCESS);
557: }

559: /* ------------------------------------------------------------------- */
560: /*
561:    ComputeInitialGuess - Calculates the initial guess

563:    Input Parameters:
564: .  user - user-defined application context
565: .  X - vector for initial guess

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

578:   PetscFunctionBeginUser;
579:   PetscCall(SNESGetDM(snes, &da));
580:   PetscCall(SNESGetApplicationContext(snes, &user));

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

585:   /* Get pointers to vector data */
586:   PetscCall(DMDAVecGetArray(da, X, &x));
587:   /* Perform local computations */
588:   for (j = ys; j < ys + ym; j++) {
589:     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;
590:   }
591:   /* Restore vectors */
592:   PetscCall(DMDAVecRestoreArray(da, X, &x));
593:   PetscFunctionReturn(PETSC_SUCCESS);
594: }

596: /*TEST

598:    test:
599:       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
600:       requires: !single

602:    test:
603:       suffix: 2
604:       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
605:       requires: !single

607: TEST*/