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 complementarity -- 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, void **);
     49: extern PetscErrorCode DestroyBoundaryConditions(void **);
     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, NULL, 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, (PetscCtxDestroyFn *)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, void **inctx)
    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, **ouser = (AppCtx **)inctx;
    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(void **ouser)
    547: {
    548:   AppCtx *user = (AppCtx *)*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:
    608:       suffix: 3
    609:       args: -snes_type vinewtonrsls -snes_linesearch_type bisection -snes_linesearch_monitor -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
    610:       requires: !single

    612:    test:
    613:       suffix: 4
    614:       args: -snes_type vinewtonssls -snes_linesearch_type bisection -snes_linesearch_monitor -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
    615:       requires: !single

    617:    test:
    618:       suffix: 5
    619:       args: -snes_type vinewtonssls -snes_linesearch_type cp -snes_linesearch_order 3 -snes_linesearch_monitor -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 -snes_linesearch_maxlambda 2.0
    620:       requires: !single

    622:    test:
    623:       suffix: 6
    624:       args: -snes_type vinewtonssls -snes_linesearch_type secant -snes_linesearch_max_it 10 -snes_linesearch_ltol 1e-2 -snes_linesearch_atol 1e-32 -snes_linesearch_monitor -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 -snes_linesearch_maxlambda 1.0
    625:       requires: !single

    627: TEST*/