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:      This is a new version of the ../tests/ex8.c code

 33:      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:      Or to run with grid sequencing on the nonlinear problem (note that you do not need to provide the number of
 36:          multigrid levels, it will be determined automatically based on the number of refinements done)

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

 41: */

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

 48: /* -------- User-defined Routines --------- */

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

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

 64:   PetscFunctionBeginUser;
 65:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));

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

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

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

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

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

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

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

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

 93:   PetscCall(SNESSetFromOptions(snes));

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

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

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

107:   PetscCall(PetscFinalize());
108:   return 0;
109: }

111: /* -------------------------------------------------------------------- */

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

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

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;

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

133: /* -------------------------------------------------------------------- */

135: /*  FormGradient - Evaluates gradient of f.

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

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;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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(PETSC_SUCCESS);
257: }

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

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

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

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;

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);
292:   hy    = 1.0 / (my + 1);
293:   hydhx = hy / hx;
294:   hxdhy = hx / hy;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

447:   PetscCall(PetscLogFlops(199.0 * mx * my));
448:   PetscFunctionReturn(PETSC_SUCCESS);
449: }

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

580:   PetscFunctionBeginUser;
581:   PetscCall(SNESGetDM(snes, &da));
582:   PetscCall(SNESGetApplicationContext(snes, &user));

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

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

598: /*TEST

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

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

609: TEST*/