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, (char *)0, 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: }
104: /* ------------------------------------------------------------------- */
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: /*
175:    RHSJacobian - User-provided routine to compute the Jacobian of
176:    the nonlinear right-hand-side function of the ODE.

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

184:    Output Parameters:
185:    J - Jacobian matrix
186:    Jpre - optionally different preconditioning matrix
187:    str - flag indicating matrix structure
188: */
189: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat J, Mat Jpre, void *ctx)
190: {
191:   DM            da;
192:   DMDALocalInfo info;
193:   PetscInt      i, j;
194:   PetscReal     hx, hy, sx, sy;

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

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

252:   PetscFunctionBeginUser;
253:   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));

255:   hx = 1.0 / (PetscReal)(Mx - 1);
256:   hy = 1.0 / (PetscReal)(My - 1);

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

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

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

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

280: /*TEST

282:     test:
283:       args: -ts_max_steps 5 -ts_monitor

285:     test:
286:       suffix: 2
287:       args: -ts_max_steps 5 -ts_monitor

289:     test:
290:       suffix: 3
291:       args: -ts_max_steps 5 -snes_fd_color -ts_monitor

293: TEST*/