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