Actual source code: ex4.c

  1: /*
  2:        The Problem:
  3:            Solve the convection-diffusion equation:

  5:              u_t+a*(u_x+u_y)=epsilon*(u_xx+u_yy)
  6:              u=0   at x=0, y=0
  7:              u_x=0 at x=1
  8:              u_y=0 at y=1
  9:              u = exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0))) at t=0

 11:        This program tests the routine of computing the Jacobian by the
 12:        finite difference method as well as PETSc.

 14: */

 16: static char help[] = "Solve the convection-diffusion equation. \n\n";

 18: #include <petscts.h>

 20: typedef struct {
 21:   PetscInt  m;       /* the number of mesh points in x-direction */
 22:   PetscInt  n;       /* the number of mesh points in y-direction */
 23:   PetscReal dx;      /* the grid space in x-direction */
 24:   PetscReal dy;      /* the grid space in y-direction */
 25:   PetscReal a;       /* the convection coefficient    */
 26:   PetscReal epsilon; /* the diffusion coefficient     */
 27:   PetscReal tfinal;
 28: } Data;

 30: extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
 31: extern PetscErrorCode Initial(Vec, void *);
 32: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
 33: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
 34: extern PetscErrorCode PostStep(TS);

 36: int main(int argc, char **argv)
 37: {
 38:   PetscInt      time_steps = 100, iout, NOUT = 1;
 39:   Vec           global;
 40:   PetscReal     dt, ftime, ftime_original;
 41:   TS            ts;
 42:   PetscViewer   viewfile;
 43:   Mat           J = 0;
 44:   Vec           x;
 45:   Data          data;
 46:   PetscInt      mn;
 47:   PetscBool     flg;
 48:   MatColoring   mc;
 49:   ISColoring    iscoloring;
 50:   MatFDColoring matfdcoloring        = 0;
 51:   PetscBool     fd_jacobian_coloring = PETSC_FALSE;
 52:   SNES          snes;
 53:   KSP           ksp;
 54:   PC            pc;

 57:   PetscInitialize(&argc, &argv, (char *)0, help);

 59:   /* set data */
 60:   data.m       = 9;
 61:   data.n       = 9;
 62:   data.a       = 1.0;
 63:   data.epsilon = 0.1;
 64:   data.dx      = 1.0 / (data.m + 1.0);
 65:   data.dy      = 1.0 / (data.n + 1.0);
 66:   mn           = (data.m) * (data.n);
 67:   PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL);

 69:   /* set initial conditions */
 70:   VecCreate(PETSC_COMM_WORLD, &global);
 71:   VecSetSizes(global, PETSC_DECIDE, mn);
 72:   VecSetFromOptions(global);
 73:   Initial(global, &data);
 74:   VecDuplicate(global, &x);

 76:   /* create timestep context */
 77:   TSCreate(PETSC_COMM_WORLD, &ts);
 78:   TSMonitorSet(ts, Monitor, &data, NULL);
 79:   TSSetType(ts, TSEULER);
 80:   dt             = 0.1;
 81:   ftime_original = data.tfinal = 1.0;

 83:   TSSetTimeStep(ts, dt);
 84:   TSSetMaxSteps(ts, time_steps);
 85:   TSSetMaxTime(ts, ftime_original);
 86:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
 87:   TSSetSolution(ts, global);

 89:   /* set user provided RHSFunction and RHSJacobian */
 90:   TSSetRHSFunction(ts, NULL, RHSFunction, &data);
 91:   MatCreate(PETSC_COMM_WORLD, &J);
 92:   MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, mn, mn);
 93:   MatSetFromOptions(J);
 94:   MatSeqAIJSetPreallocation(J, 5, NULL);
 95:   MatMPIAIJSetPreallocation(J, 5, NULL, 5, NULL);

 97:   PetscOptionsHasName(NULL, NULL, "-ts_fd", &flg);
 98:   if (!flg) {
 99:     TSSetRHSJacobian(ts, J, J, RHSJacobian, &data);
100:   } else {
101:     TSGetSNES(ts, &snes);
102:     PetscOptionsHasName(NULL, NULL, "-fd_color", &fd_jacobian_coloring);
103:     if (fd_jacobian_coloring) { /* Use finite differences with coloring */
104:       /* Get data structure of J */
105:       PetscBool pc_diagonal;
106:       PetscOptionsHasName(NULL, NULL, "-pc_diagonal", &pc_diagonal);
107:       if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */
108:         PetscInt    rstart, rend, i;
109:         PetscScalar zero = 0.0;
110:         MatGetOwnershipRange(J, &rstart, &rend);
111:         for (i = rstart; i < rend; i++) MatSetValues(J, 1, &i, 1, &i, &zero, INSERT_VALUES);
112:         MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
113:         MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
114:       } else {
115:         /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
116:         TSSetType(ts, TSBEULER);
117:         TSSetUp(ts);
118:         SNESComputeJacobianDefault(snes, x, J, J, ts);
119:       }

121:       /* create coloring context */
122:       MatColoringCreate(J, &mc);
123:       MatColoringSetType(mc, MATCOLORINGSL);
124:       MatColoringSetFromOptions(mc);
125:       MatColoringApply(mc, &iscoloring);
126:       MatColoringDestroy(&mc);
127:       MatFDColoringCreate(J, iscoloring, &matfdcoloring);
128:       MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, ts);
129:       MatFDColoringSetFromOptions(matfdcoloring);
130:       MatFDColoringSetUp(J, iscoloring, matfdcoloring);
131:       SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring);
132:       ISColoringDestroy(&iscoloring);
133:     } else { /* Use finite differences (slow) */
134:       SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, NULL);
135:     }
136:   }

138:   /* Pick up a Petsc preconditioner */
139:   /* one can always set method or preconditioner during the run time */
140:   TSGetSNES(ts, &snes);
141:   SNESGetKSP(snes, &ksp);
142:   KSPGetPC(ksp, &pc);
143:   PCSetType(pc, PCJACOBI);
144:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);

146:   TSSetFromOptions(ts);
147:   TSSetUp(ts);

149:   /* Test TSSetPostStep() */
150:   PetscOptionsHasName(NULL, NULL, "-test_PostStep", &flg);
151:   if (flg) TSSetPostStep(ts, PostStep);

153:   PetscOptionsGetInt(NULL, NULL, "-NOUT", &NOUT, NULL);
154:   for (iout = 1; iout <= NOUT; iout++) {
155:     TSSetMaxSteps(ts, time_steps);
156:     TSSetMaxTime(ts, iout * ftime_original / NOUT);
157:     TSSolve(ts, global);
158:     TSGetSolveTime(ts, &ftime);
159:     TSSetTime(ts, ftime);
160:     TSSetTimeStep(ts, dt);
161:   }
162:   /* Interpolate solution at tfinal */
163:   TSGetSolution(ts, &global);
164:   TSInterpolate(ts, ftime_original, global);

166:   PetscOptionsHasName(NULL, NULL, "-matlab_view", &flg);
167:   if (flg) { /* print solution into a MATLAB file */
168:     PetscViewerASCIIOpen(PETSC_COMM_WORLD, "out.m", &viewfile);
169:     PetscViewerPushFormat(viewfile, PETSC_VIEWER_ASCII_MATLAB);
170:     VecView(global, viewfile);
171:     PetscViewerPopFormat(viewfile);
172:     PetscViewerDestroy(&viewfile);
173:   }

175:   /* free the memories */
176:   TSDestroy(&ts);
177:   VecDestroy(&global);
178:   VecDestroy(&x);
179:   MatDestroy(&J);
180:   if (fd_jacobian_coloring) MatFDColoringDestroy(&matfdcoloring);
181:   PetscFinalize();
182:   return 0;
183: }

185: /* -------------------------------------------------------------------*/
186: /* the initial function */
187: PetscReal f_ini(PetscReal x, PetscReal y)
188: {
189:   PetscReal f;

191:   f = PetscExpReal(-20.0 * (PetscPowRealInt(x - 0.5, 2) + PetscPowRealInt(y - 0.5, 2)));
192:   return f;
193: }

195: PetscErrorCode Initial(Vec global, void *ctx)
196: {
197:   Data        *data = (Data *)ctx;
198:   PetscInt     m, row, col;
199:   PetscReal    x, y, dx, dy;
200:   PetscScalar *localptr;
201:   PetscInt     i, mybase, myend, locsize;

204:   /* make the local  copies of parameters */
205:   m  = data->m;
206:   dx = data->dx;
207:   dy = data->dy;

209:   /* determine starting point of each processor */
210:   VecGetOwnershipRange(global, &mybase, &myend);
211:   VecGetLocalSize(global, &locsize);

213:   /* Initialize the array */
214:   VecGetArrayWrite(global, &localptr);

216:   for (i = 0; i < locsize; i++) {
217:     row         = 1 + (mybase + i) - ((mybase + i) / m) * m;
218:     col         = (mybase + i) / m + 1;
219:     x           = dx * row;
220:     y           = dy * col;
221:     localptr[i] = f_ini(x, y);
222:   }

224:   VecRestoreArrayWrite(global, &localptr);
225:   return 0;
226: }

228: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, void *ctx)
229: {
230:   VecScatter         scatter;
231:   IS                 from, to;
232:   PetscInt           i, n, *idx, nsteps, maxsteps;
233:   Vec                tmp_vec;
234:   const PetscScalar *tmp;

237:   TSGetStepNumber(ts, &nsteps);
238:   /* display output at selected time steps */
239:   TSGetMaxSteps(ts, &maxsteps);
240:   if (nsteps % 10 != 0) return 0;

242:   /* Get the size of the vector */
243:   VecGetSize(global, &n);

245:   /* Set the index sets */
246:   PetscMalloc1(n, &idx);
247:   for (i = 0; i < n; i++) idx[i] = i;

249:   /* Create local sequential vectors */
250:   VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec);

252:   /* Create scatter context */
253:   ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from);
254:   ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to);
255:   VecScatterCreate(global, from, tmp_vec, to, &scatter);
256:   VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD);
257:   VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD);

259:   VecGetArrayRead(tmp_vec, &tmp);
260:   PetscPrintf(PETSC_COMM_WORLD, "At t[%" PetscInt_FMT "] =%14.2e u= %14.2e at the center \n", nsteps, (double)time, (double)PetscRealPart(tmp[n / 2]));
261:   VecRestoreArrayRead(tmp_vec, &tmp);

263:   PetscFree(idx);
264:   ISDestroy(&from);
265:   ISDestroy(&to);
266:   VecScatterDestroy(&scatter);
267:   VecDestroy(&tmp_vec);
268:   return 0;
269: }

271: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ptr)
272: {
273:   Data       *data = (Data *)ptr;
274:   PetscScalar v[5];
275:   PetscInt    idx[5], i, j, row;
276:   PetscInt    m, n, mn;
277:   PetscReal   dx, dy, a, epsilon, xc, xl, xr, yl, yr;

280:   m       = data->m;
281:   n       = data->n;
282:   mn      = m * n;
283:   dx      = data->dx;
284:   dy      = data->dy;
285:   a       = data->a;
286:   epsilon = data->epsilon;

288:   xc = -2.0 * epsilon * (1.0 / (dx * dx) + 1.0 / (dy * dy));
289:   xl = 0.5 * a / dx + epsilon / (dx * dx);
290:   xr = -0.5 * a / dx + epsilon / (dx * dx);
291:   yl = 0.5 * a / dy + epsilon / (dy * dy);
292:   yr = -0.5 * a / dy + epsilon / (dy * dy);

294:   row    = 0;
295:   v[0]   = xc;
296:   v[1]   = xr;
297:   v[2]   = yr;
298:   idx[0] = 0;
299:   idx[1] = 2;
300:   idx[2] = m;
301:   MatSetValues(A, 1, &row, 3, idx, v, INSERT_VALUES);

303:   row    = m - 1;
304:   v[0]   = 2.0 * xl;
305:   v[1]   = xc;
306:   v[2]   = yr;
307:   idx[0] = m - 2;
308:   idx[1] = m - 1;
309:   idx[2] = m - 1 + m;
310:   MatSetValues(A, 1, &row, 3, idx, v, INSERT_VALUES);

312:   for (i = 1; i < m - 1; i++) {
313:     row    = i;
314:     v[0]   = xl;
315:     v[1]   = xc;
316:     v[2]   = xr;
317:     v[3]   = yr;
318:     idx[0] = i - 1;
319:     idx[1] = i;
320:     idx[2] = i + 1;
321:     idx[3] = i + m;
322:     MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES);
323:   }

325:   for (j = 1; j < n - 1; j++) {
326:     row    = j * m;
327:     v[0]   = xc;
328:     v[1]   = xr;
329:     v[2]   = yl;
330:     v[3]   = yr;
331:     idx[0] = j * m;
332:     idx[1] = j * m;
333:     idx[2] = j * m - m;
334:     idx[3] = j * m + m;
335:     MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES);

337:     row    = j * m + m - 1;
338:     v[0]   = xc;
339:     v[1]   = 2.0 * xl;
340:     v[2]   = yl;
341:     v[3]   = yr;
342:     idx[0] = j * m + m - 1;
343:     idx[1] = j * m + m - 1 - 1;
344:     idx[2] = j * m + m - 1 - m;
345:     idx[3] = j * m + m - 1 + m;
346:     MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES);

348:     for (i = 1; i < m - 1; i++) {
349:       row    = j * m + i;
350:       v[0]   = xc;
351:       v[1]   = xl;
352:       v[2]   = xr;
353:       v[3]   = yl;
354:       v[4]   = yr;
355:       idx[0] = j * m + i;
356:       idx[1] = j * m + i - 1;
357:       idx[2] = j * m + i + 1;
358:       idx[3] = j * m + i - m;
359:       idx[4] = j * m + i + m;
360:       MatSetValues(A, 1, &row, 5, idx, v, INSERT_VALUES);
361:     }
362:   }

364:   row    = mn - m;
365:   v[0]   = xc;
366:   v[1]   = xr;
367:   v[2]   = 2.0 * yl;
368:   idx[0] = mn - m;
369:   idx[1] = mn - m + 1;
370:   idx[2] = mn - m - m;
371:   MatSetValues(A, 1, &row, 3, idx, v, INSERT_VALUES);

373:   row    = mn - 1;
374:   v[0]   = xc;
375:   v[1]   = 2.0 * xl;
376:   v[2]   = 2.0 * yl;
377:   idx[0] = mn - 1;
378:   idx[1] = mn - 2;
379:   idx[2] = mn - 1 - m;
380:   MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES);

382:   for (i = 1; i < m - 1; i++) {
383:     row    = mn - m + i;
384:     v[0]   = xl;
385:     v[1]   = xc;
386:     v[2]   = xr;
387:     v[3]   = 2.0 * yl;
388:     idx[0] = mn - m + i - 1;
389:     idx[1] = mn - m + i;
390:     idx[2] = mn - m + i + 1;
391:     idx[3] = mn - m + i - m;
392:     MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES);
393:   }

395:   MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
396:   MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);

398:   return 0;
399: }

401: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
402: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
403: {
404:   Data              *data = (Data *)ctx;
405:   PetscInt           m, n, mn;
406:   PetscReal          dx, dy;
407:   PetscReal          xc, xl, xr, yl, yr;
408:   PetscReal          a, epsilon;
409:   PetscScalar       *outptr;
410:   const PetscScalar *inptr;
411:   PetscInt           i, j, len;
412:   IS                 from, to;
413:   PetscInt          *idx;
414:   VecScatter         scatter;
415:   Vec                tmp_in, tmp_out;

418:   m       = data->m;
419:   n       = data->n;
420:   mn      = m * n;
421:   dx      = data->dx;
422:   dy      = data->dy;
423:   a       = data->a;
424:   epsilon = data->epsilon;

426:   xc = -2.0 * epsilon * (1.0 / (dx * dx) + 1.0 / (dy * dy));
427:   xl = 0.5 * a / dx + epsilon / (dx * dx);
428:   xr = -0.5 * a / dx + epsilon / (dx * dx);
429:   yl = 0.5 * a / dy + epsilon / (dy * dy);
430:   yr = -0.5 * a / dy + epsilon / (dy * dy);

432:   /* Get the length of parallel vector */
433:   VecGetSize(globalin, &len);

435:   /* Set the index sets */
436:   PetscMalloc1(len, &idx);
437:   for (i = 0; i < len; i++) idx[i] = i;

439:   /* Create local sequential vectors */
440:   VecCreateSeq(PETSC_COMM_SELF, len, &tmp_in);
441:   VecDuplicate(tmp_in, &tmp_out);

443:   /* Create scatter context */
444:   ISCreateGeneral(PETSC_COMM_SELF, len, idx, PETSC_COPY_VALUES, &from);
445:   ISCreateGeneral(PETSC_COMM_SELF, len, idx, PETSC_COPY_VALUES, &to);
446:   VecScatterCreate(globalin, from, tmp_in, to, &scatter);
447:   VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD);
448:   VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD);
449:   VecScatterDestroy(&scatter);

451:   /*Extract income array - include ghost points */
452:   VecGetArrayRead(tmp_in, &inptr);

454:   /* Extract outcome array*/
455:   VecGetArrayWrite(tmp_out, &outptr);

457:   outptr[0]     = xc * inptr[0] + xr * inptr[1] + yr * inptr[m];
458:   outptr[m - 1] = 2.0 * xl * inptr[m - 2] + xc * inptr[m - 1] + yr * inptr[m - 1 + m];
459:   for (i = 1; i < m - 1; i++) outptr[i] = xc * inptr[i] + xl * inptr[i - 1] + xr * inptr[i + 1] + yr * inptr[i + m];

461:   for (j = 1; j < n - 1; j++) {
462:     outptr[j * m]         = xc * inptr[j * m] + xr * inptr[j * m + 1] + yl * inptr[j * m - m] + yr * inptr[j * m + m];
463:     outptr[j * m + m - 1] = xc * inptr[j * m + m - 1] + 2.0 * xl * inptr[j * m + m - 1 - 1] + yl * inptr[j * m + m - 1 - m] + yr * inptr[j * m + m - 1 + m];
464:     for (i = 1; i < m - 1; i++) outptr[j * m + i] = xc * inptr[j * m + i] + xl * inptr[j * m + i - 1] + xr * inptr[j * m + i + 1] + yl * inptr[j * m + i - m] + yr * inptr[j * m + i + m];
465:   }

467:   outptr[mn - m] = xc * inptr[mn - m] + xr * inptr[mn - m + 1] + 2.0 * yl * inptr[mn - m - m];
468:   outptr[mn - 1] = 2.0 * xl * inptr[mn - 2] + xc * inptr[mn - 1] + 2.0 * yl * inptr[mn - 1 - m];
469:   for (i = 1; i < m - 1; i++) outptr[mn - m + i] = xc * inptr[mn - m + i] + xl * inptr[mn - m + i - 1] + xr * inptr[mn - m + i + 1] + 2 * yl * inptr[mn - m + i - m];

471:   VecRestoreArrayRead(tmp_in, &inptr);
472:   VecRestoreArrayWrite(tmp_out, &outptr);

474:   VecScatterCreate(tmp_out, from, globalout, to, &scatter);
475:   VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD);
476:   VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD);

478:   /* Destroy idx aand scatter */
479:   VecDestroy(&tmp_in);
480:   VecDestroy(&tmp_out);
481:   ISDestroy(&from);
482:   ISDestroy(&to);
483:   VecScatterDestroy(&scatter);

485:   PetscFree(idx);
486:   return 0;
487: }

489: PetscErrorCode PostStep(TS ts)
490: {
491:   PetscReal t;

494:   TSGetTime(ts, &t);
495:   PetscPrintf(PETSC_COMM_SELF, "  PostStep, t: %g\n", (double)t);
496:   return 0;
497: }

499: /*TEST

501:     test:
502:       args: -ts_fd -ts_type beuler
503:       output_file: output/ex4.out

505:     test:
506:       suffix: 2
507:       args: -ts_fd -ts_type beuler
508:       nsize: 2
509:       output_file: output/ex4.out

511:     test:
512:       suffix: 3
513:       args: -ts_fd -ts_type cn

515:     test:
516:       suffix: 4
517:       args: -ts_fd -ts_type cn
518:       output_file: output/ex4_3.out
519:       nsize: 2

521:     test:
522:       suffix: 5
523:       args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
524:       output_file: output/ex4.out

526:     test:
527:       suffix: 6
528:       args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
529:       output_file: output/ex4.out
530:       nsize: 2

532:     test:
533:       suffix: 7
534:       requires: !single
535:       args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1

537:     test:
538:       suffix: 8
539:       requires: !single
540:       args: -ts_type rk -ts_rk_type 5dp -ts_dt .01 -ts_adapt_type none -ts_view

542: TEST*/