Actual source code: ex13.c

      1: static char help[] = "Time-dependent PDE in 2d. Simplified from ex7.c for illustrating how to use TS on a structured domain. \n";
      2: /*
      3:    u_t = uxx + uyy
      4:    0 < x < 1, 0 < y < 1;
      5:    At t=0: u(x,y) = exp(c*r*r*r), if r=PetscSqrtReal((x-.5)*(x-.5) + (y-.5)*(y-.5)) < .125
      6:            u(x,y) = 0.0           if r >= .125

      8:     mpiexec -n 2 ./ex13 -da_grid_x 40 -da_grid_y 40 -ts_max_steps 2 -snes_monitor -ksp_monitor
      9:     mpiexec -n 1 ./ex13 -snes_fd_color -ts_monitor_draw_solution
     10:     mpiexec -n 2 ./ex13 -ts_type sundials -ts_monitor
     11: */

     13: #include <petscdm.h>
     14: #include <petscdmda.h>
     15: #include <petscts.h>

     17: /*
     18:    User-defined data structures and routines
     19: */
     20: typedef struct {
     21:   PetscReal c;
     22: } AppCtx;

     24: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
     25: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
     26: extern PetscErrorCode FormInitialSolution(DM, Vec, void *);

     28: int main(int argc, char **argv)
     29: {
     30:   TS        ts;    /* nonlinear solver */
     31:   Vec       u, r;  /* solution, residual vector */
     32:   Mat       J;     /* Jacobian matrix */
     33:   PetscInt  steps; /* iterations for convergence */
     34:   DM        da;
     35:   PetscReal ftime, dt;
     36:   AppCtx    user; /* user-defined work context */

     38:   PetscFunctionBeginUser;
     39:   PetscCall(PetscInitialize(&argc, &argv, NULL, help));
     40:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     41:      Create distributed array (DMDA) to manage parallel grid and vectors
     42:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     43:   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));
     44:   PetscCall(DMSetFromOptions(da));
     45:   PetscCall(DMSetUp(da));

     47:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     48:      Extract global vectors from DMDA;
     49:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     50:   PetscCall(DMCreateGlobalVector(da, &u));
     51:   PetscCall(VecDuplicate(u, &r));

     53:   /* Initialize user application context */
     54:   user.c = -30.0;

     56:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     57:      Create timestepping solver context
     58:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     59:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
     60:   PetscCall(TSSetDM(ts, da));
     61:   PetscCall(TSSetType(ts, TSBEULER));
     62:   PetscCall(TSSetRHSFunction(ts, r, RHSFunction, &user));

     64:   /* Set Jacobian */
     65:   PetscCall(DMSetMatType(da, MATAIJ));
     66:   PetscCall(DMCreateMatrix(da, &J));
     67:   PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, NULL));

     69:   ftime = 1.0;
     70:   PetscCall(TSSetMaxTime(ts, ftime));
     71:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));

     73:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     74:      Set initial conditions
     75:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     76:   PetscCall(FormInitialSolution(da, u, &user));
     77:   dt = .01;
     78:   PetscCall(TSSetTimeStep(ts, dt));

     80:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     81:      Set runtime options
     82:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     83:   PetscCall(TSSetFromOptions(ts));

     85:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     86:      Solve nonlinear system
     87:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     88:   PetscCall(TSSolve(ts, u));
     89:   PetscCall(TSGetSolveTime(ts, &ftime));
     90:   PetscCall(TSGetStepNumber(ts, &steps));

     92:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
     93:      Free work space.
     94:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
     95:   PetscCall(MatDestroy(&J));
     96:   PetscCall(VecDestroy(&u));
     97:   PetscCall(VecDestroy(&r));
     98:   PetscCall(TSDestroy(&ts));
     99:   PetscCall(DMDestroy(&da));

    101:   PetscCall(PetscFinalize());
    102:   return 0;
    103: }

    105: /*
    106:    RHSFunction - Evaluates nonlinear function, F(u).

    108:    Input Parameters:
    109: .  ts - the TS context
    110: .  U - input vector
    111: .  ptr - optional user-defined context, as set by TSSetFunction()

    113:    Output Parameter:
    114: .  F - function vector
    115:  */
    116: PetscErrorCode RHSFunction(TS ts, PetscReal ftime, Vec U, Vec F, void *ptr)
    117: {
    118:   /* PETSC_UNUSED AppCtx *user=(AppCtx*)ptr; */
    119:   DM          da;
    120:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
    121:   PetscReal   two = 2.0, hx, hy, sx, sy;
    122:   PetscScalar u, uxx, uyy, **uarray, **f;
    123:   Vec         localU;

    125:   PetscFunctionBeginUser;
    126:   PetscCall(TSGetDM(ts, &da));
    127:   PetscCall(DMGetLocalVector(da, &localU));
    128:   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));

    130:   hx = 1.0 / (PetscReal)(Mx - 1);
    131:   sx = 1.0 / (hx * hx);
    132:   hy = 1.0 / (PetscReal)(My - 1);
    133:   sy = 1.0 / (hy * hy);

    135:   /*
    136:      Scatter ghost points to local vector,using the 2-step process
    137:         DMGlobalToLocalBegin(),DMGlobalToLocalEnd().
    138:      By placing code between these two statements, computations can be
    139:      done while messages are in transition.
    140:   */
    141:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
    142:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));

    144:   /* Get pointers to vector data */
    145:   PetscCall(DMDAVecGetArrayRead(da, localU, &uarray));
    146:   PetscCall(DMDAVecGetArray(da, F, &f));

    148:   /* Get local grid boundaries */
    149:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

    151:   /* Compute function over the locally owned part of the grid */
    152:   for (j = ys; j < ys + ym; j++) {
    153:     for (i = xs; i < xs + xm; i++) {
    154:       if (i == 0 || j == 0 || i == Mx - 1 || j == My - 1) {
    155:         f[j][i] = uarray[j][i];
    156:         continue;
    157:       }
    158:       u       = uarray[j][i];
    159:       uxx     = (-two * u + uarray[j][i - 1] + uarray[j][i + 1]) * sx;
    160:       uyy     = (-two * u + uarray[j - 1][i] + uarray[j + 1][i]) * sy;
    161:       f[j][i] = uxx + uyy;
    162:     }
    163:   }

    165:   /* Restore vectors */
    166:   PetscCall(DMDAVecRestoreArrayRead(da, localU, &uarray));
    167:   PetscCall(DMDAVecRestoreArray(da, F, &f));
    168:   PetscCall(DMRestoreLocalVector(da, &localU));
    169:   PetscCall(PetscLogFlops(11.0 * ym * xm));
    170:   PetscFunctionReturn(PETSC_SUCCESS);
    171: }

    173: /*
    174:    RHSJacobian - User-provided routine to compute the Jacobian of
    175:    the nonlinear right-hand-side function of the ODE.

    177:    Input Parameters:
    178:    ts - the TS context
    179:    t - current time
    180:    U - global input vector
    181:    dummy - optional user-defined context, as set by TSetRHSJacobian()

    183:    Output Parameters:
    184:    J - Jacobian matrix
    185:    Jpre - optionally different matrix used to construct the preconditioner

    187: */
    188: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx)
    189: {
    190:   DM            da;
    191:   DMDALocalInfo info;
    192:   PetscInt      i, j;
    193:   PetscReal     hx, hy, sx, sy;

    195:   PetscFunctionBeginUser;
    196:   PetscCall(TSGetDM(ts, &da));
    197:   PetscCall(DMDAGetLocalInfo(da, &info));
    198:   hx = 1.0 / (PetscReal)(info.mx - 1);
    199:   sx = 1.0 / (hx * hx);
    200:   hy = 1.0 / (PetscReal)(info.my - 1);
    201:   sy = 1.0 / (hy * hy);
    202:   for (j = info.ys; j < info.ys + info.ym; j++) {
    203:     for (i = info.xs; i < info.xs + info.xm; i++) {
    204:       PetscInt    nc = 0;
    205:       MatStencil  row, col[5];
    206:       PetscScalar val[5];
    207:       row.i = i;
    208:       row.j = j;
    209:       if (i == 0 || j == 0 || i == info.mx - 1 || j == info.my - 1) {
    210:         col[nc].i = i;
    211:         col[nc].j = j;
    212:         val[nc++] = 1.0;
    213:       } else {
    214:         col[nc].i = i - 1;
    215:         col[nc].j = j;
    216:         val[nc++] = sx;
    217:         col[nc].i = i + 1;
    218:         col[nc].j = j;
    219:         val[nc++] = sx;
    220:         col[nc].i = i;
    221:         col[nc].j = j - 1;
    222:         val[nc++] = sy;
    223:         col[nc].i = i;
    224:         col[nc].j = j + 1;
    225:         val[nc++] = sy;
    226:         col[nc].i = i;
    227:         col[nc].j = j;
    228:         val[nc++] = -2 * sx - 2 * sy;
    229:       }
    230:       PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
    231:     }
    232:   }
    233:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
    234:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
    235:   if (J != Jpre) {
    236:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
    237:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
    238:   }
    239:   PetscFunctionReturn(PETSC_SUCCESS);
    240: }

    242: PetscErrorCode FormInitialSolution(DM da, Vec U, void *ptr)
    243: {
    244:   AppCtx       *user = (AppCtx *)ptr;
    245:   PetscReal     c    = user->c;
    246:   PetscInt      i, j, xs, ys, xm, ym, Mx, My;
    247:   PetscScalar **u;
    248:   PetscReal     hx, hy, x, y, r;

    250:   PetscFunctionBeginUser;
    251:   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));

    253:   hx = 1.0 / (PetscReal)(Mx - 1);
    254:   hy = 1.0 / (PetscReal)(My - 1);

    256:   /* Get pointers to vector data */
    257:   PetscCall(DMDAVecGetArray(da, U, &u));

    259:   /* Get local grid boundaries */
    260:   PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL));

    262:   /* Compute function over the locally owned part of the grid */
    263:   for (j = ys; j < ys + ym; j++) {
    264:     y = j * hy;
    265:     for (i = xs; i < xs + xm; i++) {
    266:       x = i * hx;
    267:       r = PetscSqrtReal((x - .5) * (x - .5) + (y - .5) * (y - .5));
    268:       if (r < .125) u[j][i] = PetscExpReal(c * r * r * r);
    269:       else u[j][i] = 0.0;
    270:     }
    271:   }

    273:   /* Restore vectors */
    274:   PetscCall(DMDAVecRestoreArray(da, U, &u));
    275:   PetscFunctionReturn(PETSC_SUCCESS);
    276: }

    278: /*TEST

    280:     test:
    281:       args: -ts_max_steps 5 -ts_monitor

    283:     test:
    284:       suffix: 2
    285:       args: -ts_max_steps 5 -ts_monitor

    287:     test:
    288:       suffix: 3
    289:       args: -ts_max_steps 5 -snes_fd_color -ts_monitor

    291: TEST*/