Actual source code: ex13.c


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

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

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

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

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

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

 40:   PetscInitialize(&argc, &argv, (char *)0, help);
 41:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 42:      Create distributed array (DMDA) to manage parallel grid and vectors
 43:   - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 44:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, DM_BOUNDARY_NONE, DMDA_STENCIL_STAR, 8, 8, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, NULL, &da);
 45:   DMSetFromOptions(da);
 46:   DMSetUp(da);

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

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

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

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

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

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

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

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

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

102:   PetscFinalize();
103:   return 0;
104: }
105: /* ------------------------------------------------------------------- */
106: /*
107:    RHSFunction - Evaluates nonlinear function, F(u).

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

276:   /* Restore vectors */
277:   DMDAVecRestoreArray(da, U, &u);
278:   return 0;
279: }

281: /*TEST

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

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

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

294: TEST*/