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 */

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

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

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

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

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

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

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

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

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

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

102:   PetscCall(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;

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

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:   PetscCall(DMGlobalToLocalBegin(da, U, INSERT_VALUES, localU));
143:   PetscCall(DMGlobalToLocalEnd(da, U, INSERT_VALUES, localU));

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

149:   /* Get local grid boundaries */
150:   PetscCall(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:   PetscCall(DMDAVecRestoreArrayRead(da, localU, &uarray));
168:   PetscCall(DMDAVecRestoreArray(da, F, &f));
169:   PetscCall(DMRestoreLocalVector(da, &localU));
170:   PetscCall(PetscLogFlops(11.0 * ym * xm));
171:   PetscFunctionReturn(PETSC_SUCCESS);
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;

197:   PetscFunctionBeginUser;
198:   PetscCall(TSGetDM(ts, &da));
199:   PetscCall(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:       PetscCall(MatSetValuesStencil(Jpre, 1, &row, nc, col, val, INSERT_VALUES));
233:     }
234:   }
235:   PetscCall(MatAssemblyBegin(Jpre, MAT_FINAL_ASSEMBLY));
236:   PetscCall(MatAssemblyEnd(Jpre, MAT_FINAL_ASSEMBLY));
237:   if (J != Jpre) {
238:     PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
239:     PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
240:   }
241:   PetscFunctionReturn(PETSC_SUCCESS);
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;

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

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

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

262:   /* Get local grid boundaries */
263:   PetscCall(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:   PetscCall(DMDAVecRestoreArray(da, U, &u));
278:   PetscFunctionReturn(PETSC_SUCCESS);
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*/