Actual source code: ex7.c

      1: static char help[] = "Nonlinear, time-dependent PDE in 2d.\n";

      3: /*
      4:    Include "petscdmda.h" so that we can use distributed arrays (DMDAs).
      5:    Include "petscts.h" so that we can use SNES solvers.  Note that this
      6:    file automatically includes:
      7:      petscsys.h       - base PETSc routines   petscvec.h - vectors
      8:      petscmat.h - matrices
      9:      petscis.h     - index sets            petscksp.h - Krylov subspace methods
     10:      petscviewer.h - viewers               petscpc.h  - preconditioners
     11:      petscksp.h   - linear solvers
     12: */
     13: #include <petscdm.h>
     14: #include <petscdmda.h>
     15: #include <petscts.h>

     17: /*
     18:    User-defined routines
     19: */
     20: extern PetscErrorCode FormFunction(TS, PetscReal, Vec, Vec, void *), FormInitialSolution(DM, Vec);
     21: extern PetscErrorCode MyTSMonitor(TS, PetscInt, PetscReal, Vec, void *);
     22: extern PetscErrorCode MySNESMonitor(SNES, PetscInt, PetscReal, PetscViewerAndFormat *);

     24: int main(int argc, char **argv)
     25: {
     26:   TS                    ts; /* time integrator */
     27:   SNES                  snes;
     28:   Vec                   x, r; /* solution, residual vectors */
     29:   DM                    da;
     30:   PetscViewerAndFormat *vf;

     32:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     33:      Initialize program
     34:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     35:   PetscFunctionBeginUser;
     36:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     37:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     38:      Create distributed array (DMDA) to manage parallel grid and vectors
     39:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     40:   PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 8, 8, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da));
     41:   PetscCall(DMSetFromOptions(da));
     42:   PetscCall(DMSetUp(da));

     44:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     45:      Extract global vectors from DMDA; then duplicate for remaining
     46:      vectors that are the same types
     47:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     48:   PetscCall(DMCreateGlobalVector(da, &x));
     49:   PetscCall(VecDuplicate(x, &r));

     51:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     52:      Create timestepping solver context
     53:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     54:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     55:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
     56:   PetscCall(TSSetRHSFunction(ts, NULL, FormFunction, da));

     58:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     59:      Create matrix data structure; set Jacobian evaluation routine

     61:      Set Jacobian matrix data structure and default Jacobian evaluation
     62:      routine. User can override with:
     63:      -snes_mf : matrix-free Newton-Krylov method with no preconditioning
     64:                 (unless user explicitly sets preconditioner)
     65:      -snes_mf_operator : form matrix used to construct the preconditioner as set by the user,
     66:                          but use matrix-free approx for Jacobian-vector
     67:                          products within Newton-Krylov method

     69:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */

     71:   PetscCall(TSSetMaxTime(ts, 1.0));
     72:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
     73:   PetscCall(TSMonitorSet(ts, MyTSMonitor, PETSC_VIEWER_STDOUT_WORLD, NULL));
     74:   PetscCall(TSSetDM(ts, da));
     75:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     76:      Customize nonlinear solver
     77:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     78:   PetscCall(TSSetType(ts, TSBEULER));
     79:   PetscCall(TSGetSNES(ts, &snes));
     80:   PetscCall(PetscViewerAndFormatCreate(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_DEFAULT, &vf));
     81:   PetscCall(SNESMonitorSet(snes, (PetscErrorCode (*)(SNES, PetscInt, PetscReal, void *))MySNESMonitor, vf, (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy));

     83:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     84:      Set initial conditions
     85:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     86:   PetscCall(FormInitialSolution(da, x));
     87:   PetscCall(TSSetTimeStep(ts, .0001));
     88:   PetscCall(TSSetSolution(ts, x));

     90:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     91:      Set runtime options
     92:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     93:   PetscCall(TSSetFromOptions(ts));

     95:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     96:      Solve nonlinear system
     97:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     98:   PetscCall(TSSolve(ts, x));

    100:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    101:      Free work space.  All PETSc objects should be destroyed when they
    102:      are no longer needed.
    103:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
    104:   PetscCall(VecDestroy(&x));
    105:   PetscCall(VecDestroy(&r));
    106:   PetscCall(TSDestroy(&ts));
    107:   PetscCall(DMDestroy(&da));

    109:   PetscCall(PetscFinalize());
    110:   return 0;
    111: }
    112: /* ------------------------------------------------------------------- */
    113: /*
    114:    FormFunction - Evaluates nonlinear function, F(x).

    116:    Input Parameters:
    117: .  ts - the TS context
    118: .  X - input vector
    119: .  ptr - optional user-defined context, as set by SNESSetFunction()

    121:    Output Parameter:
    122: .  F - function vector
    123:  */
    124: PetscErrorCode FormFunction(TS ts, PetscReal ftime, Vec X, Vec F, void *ptr)
    125: {
    126:   DM          da;
    127:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
    128:   PetscReal   two = 2.0, hx, hy, sx, sy;
    129:   PetscScalar u, uxx, uyy, **x, **f;
    130:   Vec         localX;

    132:   PetscFunctionBeginUser;
    133:   PetscCall(TSGetDM(ts, &da));
    134:   PetscCall(DMGetLocalVector(da, &localX));
    135:   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));

    137:   hx = 1.0 / (PetscReal)(Mx - 1);
    138:   sx = 1.0 / (hx * hx);
    139:   hy = 1.0 / (PetscReal)(My - 1);
    140:   sy = 1.0 / (hy * hy);

    142:   /*
    143:      Scatter ghost points to local vector,using the 2-step process
    144:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    145:      By placing code between these two statements, computations can be
    146:      done while messages are in transition.
    147:   */
    148:   PetscCall(DMGlobalToLocalBegin(da, X, INSERT_VALUES, localX));
    149:   PetscCall(DMGlobalToLocalEnd(da, X, INSERT_VALUES, localX));

    151:   /*
    152:      Get pointers to vector data
    153:   */
    154:   PetscCall(DMDAVecGetArrayRead(da, localX, &x));
    155:   PetscCall(DMDAVecGetArray(da, F, &f));

    157:   /*
    158:      Get local grid boundaries
    159:   */
    160:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

    162:   /*
    163:      Compute function over the locally owned part of the grid
    164:   */
    165:   for (j = ys; j < ys + ym; j++) {
    166:     for (i = xs; i < xs + xm; i++) {
    167:       if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
    168:         f[j][i] = x[j][i];
    169:         continue;
    170:       }
    171:       u   = x[j][i];
    172:       uxx = (two * u - x[j][i - 1] - x[j][i + 1]) * sx;
    173:       uyy = (two * u - x[j - 1][i] - x[j + 1][i]) * sy;
    174:       /*      f[j][i] = -(uxx + uyy); */
    175:       f[j][i] = -u * (uxx + uyy) - (4.0 - 1.0) * ((x[j][i + 1] - x[j][i - 1]) * (x[j][i + 1] - x[j][i - 1]) * .25 * sx + (x[j + 1][i] - x[j - 1][i]) * (x[j + 1][i] - x[j - 1][i]) * .25 * sy);
    176:     }
    177:   }

    179:   /*
    180:      Restore vectors
    181:   */
    182:   PetscCall(DMDAVecRestoreArrayRead(da, localX, &x));
    183:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    184:   PetscCall(DMRestoreLocalVector(da, &localX));
    185:   PetscCall(PetscLogFlops(11.0 * ym * xm));
    186:   PetscFunctionReturn(PETSC_SUCCESS);
    187: }

    189: /* ------------------------------------------------------------------- */
    190: PetscErrorCode FormInitialSolution(DM da, Vec U)
    191: {
    192:   PetscInt      i, j, xs, ys, xm, ym, Mx, My;
    193:   PetscScalar **u;
    194:   PetscReal     hx, hy, x, y, r;

    196:   PetscFunctionBeginUser;
    197:   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));

    199:   hx = 1.0 / (PetscReal)(Mx - 1);
    200:   hy = 1.0 / (PetscReal)(My - 1);

    202:   /*
    203:      Get pointers to vector data
    204:   */
    205:   PetscCall(DMDAVecGetArray(da, U, &u));

    207:   /*
    208:      Get local grid boundaries
    209:   */
    210:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

    212:   /*
    213:      Compute function over the locally owned part of the grid
    214:   */
    215:   for (j = ys; j < ys + ym; j++) {
    216:     y = j * hy;
    217:     for (i = xs; i < xs + xm; i++) {
    218:       x = i * hx;
    219:       r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
    220:       if (r < .125) u[j][i] = PetscExpReal(-30.0 * r * r * r);
    221:       else u[j][i] = 0.0;
    222:     }
    223:   }

    225:   /*
    226:      Restore vectors
    227:   */
    228:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    229:   PetscFunctionReturn(PETSC_SUCCESS);
    230: }

    232: PetscErrorCode MyTSMonitor(TS ts, PetscInt step, PetscReal ptime, Vec v, void *ctx)
    233: {
    234:   PetscReal norm;
    235:   MPI_Comm  comm;

    237:   PetscFunctionBeginUser;
    238:   if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* step of -1 indicates an interpolated solution */
    239:   PetscCall(VecNorm(v, NORM_2, &norm));
    240:   PetscCall(PetscObjectGetComm((PetscObject)ts, &comm));
    241:   PetscCall(PetscPrintf(comm, "timestep %" PetscInt_FMT " time %g norm %g\n", step, (double)ptime, (double)norm));
    242:   PetscFunctionReturn(PETSC_SUCCESS);
    243: }

    245: /*
    246:    MySNESMonitor - illustrate how to set user-defined monitoring routine for SNES.
    247:    Input Parameters:
    248:      snes - the SNES context
    249:      its - iteration number
    250:      fnorm - 2-norm function value (may be estimated)
    251:      ctx - optional user-defined context for private data for the
    252:          monitor routine, as set by SNESMonitorSet()
    253:  */
    254: PetscErrorCode MySNESMonitor(SNES snes, PetscInt its, PetscReal fnorm, PetscViewerAndFormat *vf)
    255: {
    256:   PetscFunctionBeginUser;
    257:   PetscCall(SNESMonitorDefaultShort(snes, its, fnorm, vf));
    258:   PetscFunctionReturn(PETSC_SUCCESS);
    259: }

    261: /*TEST

    263:     test:
    264:       args: -ts_max_steps 5

    266:     test:
    267:       suffix: 2
    268:       args: -ts_max_steps 5 -snes_mf_operator

    270:     test:
    271:       suffix: 3
    272:       args: -ts_max_steps 5 -snes_mf -pc_type none

    274: TEST*/