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;
56: PetscFunctionBeginUser;
57: PetscCall(PetscInitialize(&argc, &argv, NULL, 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: PetscCall(PetscOptionsGetInt(NULL, NULL, "-time", &time_steps, NULL));
69: /* set initial conditions */
70: PetscCall(VecCreate(PETSC_COMM_WORLD, &global));
71: PetscCall(VecSetSizes(global, PETSC_DECIDE, mn));
72: PetscCall(VecSetFromOptions(global));
73: PetscCall(Initial(global, &data));
74: PetscCall(VecDuplicate(global, &x));
76: /* create timestep context */
77: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
78: PetscCall(TSMonitorSet(ts, Monitor, &data, NULL));
79: PetscCall(TSSetType(ts, TSEULER));
80: dt = 0.1;
81: ftime_original = data.tfinal = 1.0;
83: PetscCall(TSSetTimeStep(ts, dt));
84: PetscCall(TSSetMaxSteps(ts, time_steps));
85: PetscCall(TSSetMaxTime(ts, ftime_original));
86: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
87: PetscCall(TSSetSolution(ts, global));
89: /* set user provided RHSFunction and RHSJacobian */
90: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &data));
91: PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
92: PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, mn, mn));
93: PetscCall(MatSetFromOptions(J));
94: PetscCall(MatSeqAIJSetPreallocation(J, 5, NULL));
95: PetscCall(MatMPIAIJSetPreallocation(J, 5, NULL, 5, NULL));
97: PetscCall(PetscOptionsHasName(NULL, NULL, "-ts_fd", &flg));
98: if (!flg) PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, &data));
99: else {
100: PetscCall(TSGetSNES(ts, &snes));
101: PetscCall(PetscOptionsHasName(NULL, NULL, "-fd_color", &fd_jacobian_coloring));
102: if (fd_jacobian_coloring) { /* Use finite differences with coloring */
103: /* Get data structure of J */
104: PetscBool pc_diagonal;
105: PetscCall(PetscOptionsHasName(NULL, NULL, "-pc_diagonal", &pc_diagonal));
106: if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */
107: PetscInt rstart, rend, i;
108: PetscScalar zero = 0.0;
109: PetscCall(MatGetOwnershipRange(J, &rstart, &rend));
110: for (i = rstart; i < rend; i++) PetscCall(MatSetValues(J, 1, &i, 1, &i, &zero, INSERT_VALUES));
111: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
112: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
113: } else {
114: /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
115: PetscCall(TSSetType(ts, TSBEULER));
116: PetscCall(TSSetUp(ts));
117: PetscCall(SNESComputeJacobianDefault(snes, x, J, J, ts));
118: }
120: /* create coloring context */
121: PetscCall(MatColoringCreate(J, &mc));
122: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
123: PetscCall(MatColoringSetFromOptions(mc));
124: PetscCall(MatColoringApply(mc, &iscoloring));
125: PetscCall(MatColoringDestroy(&mc));
126: PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring));
127: PetscCall(MatFDColoringSetFunction(matfdcoloring, (MatFDColoringFn *)SNESTSFormFunction, ts));
128: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
129: PetscCall(MatFDColoringSetUp(J, iscoloring, matfdcoloring));
130: PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring));
131: PetscCall(ISColoringDestroy(&iscoloring));
132: } else { /* Use finite differences (slow) */
133: PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, NULL));
134: }
135: }
137: /* Pick up a PETSc preconditioner */
138: /* one can always set method or preconditioner during the run time */
139: PetscCall(TSGetSNES(ts, &snes));
140: PetscCall(SNESGetKSP(snes, &ksp));
141: PetscCall(KSPGetPC(ksp, &pc));
142: PetscCall(PCSetType(pc, PCJACOBI));
143: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
145: PetscCall(TSSetFromOptions(ts));
146: PetscCall(TSSetUp(ts));
148: /* Test TSSetPostStep() */
149: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_PostStep", &flg));
150: if (flg) PetscCall(TSSetPostStep(ts, PostStep));
152: PetscCall(PetscOptionsGetInt(NULL, NULL, "-NOUT", &NOUT, NULL));
153: for (iout = 1; iout <= NOUT; iout++) {
154: PetscCall(TSSetMaxSteps(ts, time_steps));
155: PetscCall(TSSetMaxTime(ts, iout * ftime_original / NOUT));
156: PetscCall(TSSolve(ts, global));
157: PetscCall(TSGetSolveTime(ts, &ftime));
158: PetscCall(TSSetTime(ts, ftime));
159: PetscCall(TSSetTimeStep(ts, dt));
160: }
161: /* Interpolate solution at tfinal */
162: PetscCall(TSGetSolution(ts, &global));
163: PetscCall(TSInterpolate(ts, ftime_original, global));
165: PetscCall(PetscOptionsHasName(NULL, NULL, "-matlab_view", &flg));
166: if (flg) { /* print solution into a MATLAB file */
167: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "out.m", &viewfile));
168: PetscCall(PetscViewerPushFormat(viewfile, PETSC_VIEWER_ASCII_MATLAB));
169: PetscCall(VecView(global, viewfile));
170: PetscCall(PetscViewerPopFormat(viewfile));
171: PetscCall(PetscViewerDestroy(&viewfile));
172: }
174: /* free the memories */
175: PetscCall(TSDestroy(&ts));
176: PetscCall(VecDestroy(&global));
177: PetscCall(VecDestroy(&x));
178: PetscCall(MatDestroy(&J));
179: if (fd_jacobian_coloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
180: PetscCall(PetscFinalize());
181: return 0;
182: }
184: /* -------------------------------------------------------------------*/
185: /* the initial function */
186: PetscReal f_ini(PetscReal x, PetscReal y)
187: {
188: PetscReal f;
190: f = PetscExpReal(-20.0 * (PetscPowRealInt(x - 0.5, 2) + PetscPowRealInt(y - 0.5, 2)));
191: return f;
192: }
194: PetscErrorCode Initial(Vec global, PetscCtx ctx)
195: {
196: Data *data = (Data *)ctx;
197: PetscInt m, row, col;
198: PetscReal x, y, dx, dy;
199: PetscScalar *localptr;
200: PetscInt i, mybase, myend, locsize;
202: PetscFunctionBeginUser;
203: /* make the local copies of parameters */
204: m = data->m;
205: dx = data->dx;
206: dy = data->dy;
208: /* determine starting point of each processor */
209: PetscCall(VecGetOwnershipRange(global, &mybase, &myend));
210: PetscCall(VecGetLocalSize(global, &locsize));
212: /* Initialize the array */
213: PetscCall(VecGetArrayWrite(global, &localptr));
215: for (i = 0; i < locsize; i++) {
216: row = 1 + (mybase + i) - ((mybase + i) / m) * m;
217: col = (mybase + i) / m + 1;
218: x = dx * row;
219: y = dy * col;
220: localptr[i] = f_ini(x, y);
221: }
223: PetscCall(VecRestoreArrayWrite(global, &localptr));
224: PetscFunctionReturn(PETSC_SUCCESS);
225: }
227: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec global, PetscCtx ctx)
228: {
229: VecScatter scatter;
230: IS from, to;
231: PetscInt i, n, *idx, nsteps, maxsteps;
232: Vec tmp_vec;
233: const PetscScalar *tmp;
235: PetscFunctionBeginUser;
236: PetscCall(TSGetStepNumber(ts, &nsteps));
237: /* display output at selected time steps */
238: PetscCall(TSGetMaxSteps(ts, &maxsteps));
239: if (nsteps % 10 != 0) PetscFunctionReturn(PETSC_SUCCESS);
241: /* Get the size of the vector */
242: PetscCall(VecGetSize(global, &n));
244: /* Set the index sets */
245: PetscCall(PetscMalloc1(n, &idx));
246: for (i = 0; i < n; i++) idx[i] = i;
248: /* Create local sequential vectors */
249: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec));
251: /* Create scatter context */
252: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from));
253: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to));
254: PetscCall(VecScatterCreate(global, from, tmp_vec, to, &scatter));
255: PetscCall(VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD));
256: PetscCall(VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD));
258: PetscCall(VecGetArrayRead(tmp_vec, &tmp));
259: PetscCall(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])));
260: PetscCall(VecRestoreArrayRead(tmp_vec, &tmp));
262: PetscCall(PetscFree(idx));
263: PetscCall(ISDestroy(&from));
264: PetscCall(ISDestroy(&to));
265: PetscCall(VecScatterDestroy(&scatter));
266: PetscCall(VecDestroy(&tmp_vec));
267: PetscFunctionReturn(PETSC_SUCCESS);
268: }
270: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat BB, void *ptr)
271: {
272: Data *data = (Data *)ptr;
273: PetscScalar v[5];
274: PetscInt idx[5], i, j, row;
275: PetscInt m, n, mn;
276: PetscReal dx, dy, a, epsilon, xc, xl, xr, yl, yr;
278: PetscFunctionBeginUser;
279: m = data->m;
280: n = data->n;
281: mn = m * n;
282: dx = data->dx;
283: dy = data->dy;
284: a = data->a;
285: epsilon = data->epsilon;
287: xc = -2.0 * epsilon * (1.0 / (dx * dx) + 1.0 / (dy * dy));
288: xl = 0.5 * a / dx + epsilon / (dx * dx);
289: xr = -0.5 * a / dx + epsilon / (dx * dx);
290: yl = 0.5 * a / dy + epsilon / (dy * dy);
291: yr = -0.5 * a / dy + epsilon / (dy * dy);
293: row = 0;
294: v[0] = xc;
295: v[1] = xr;
296: v[2] = yr;
297: idx[0] = 0;
298: idx[1] = 2;
299: idx[2] = m;
300: PetscCall(MatSetValues(A, 1, &row, 3, idx, v, INSERT_VALUES));
302: row = m - 1;
303: v[0] = 2.0 * xl;
304: v[1] = xc;
305: v[2] = yr;
306: idx[0] = m - 2;
307: idx[1] = m - 1;
308: idx[2] = m - 1 + m;
309: PetscCall(MatSetValues(A, 1, &row, 3, idx, v, INSERT_VALUES));
311: for (i = 1; i < m - 1; i++) {
312: row = i;
313: v[0] = xl;
314: v[1] = xc;
315: v[2] = xr;
316: v[3] = yr;
317: idx[0] = i - 1;
318: idx[1] = i;
319: idx[2] = i + 1;
320: idx[3] = i + m;
321: PetscCall(MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES));
322: }
324: for (j = 1; j < n - 1; j++) {
325: row = j * m;
326: v[0] = xc;
327: v[1] = xr;
328: v[2] = yl;
329: v[3] = yr;
330: idx[0] = j * m;
331: idx[1] = j * m;
332: idx[2] = j * m - m;
333: idx[3] = j * m + m;
334: PetscCall(MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES));
336: row = j * m + m - 1;
337: v[0] = xc;
338: v[1] = 2.0 * xl;
339: v[2] = yl;
340: v[3] = yr;
341: idx[0] = j * m + m - 1;
342: idx[1] = j * m + m - 1 - 1;
343: idx[2] = j * m + m - 1 - m;
344: idx[3] = j * m + m - 1 + m;
345: PetscCall(MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES));
347: for (i = 1; i < m - 1; i++) {
348: row = j * m + i;
349: v[0] = xc;
350: v[1] = xl;
351: v[2] = xr;
352: v[3] = yl;
353: v[4] = yr;
354: idx[0] = j * m + i;
355: idx[1] = j * m + i - 1;
356: idx[2] = j * m + i + 1;
357: idx[3] = j * m + i - m;
358: idx[4] = j * m + i + m;
359: PetscCall(MatSetValues(A, 1, &row, 5, idx, v, INSERT_VALUES));
360: }
361: }
363: row = mn - m;
364: v[0] = xc;
365: v[1] = xr;
366: v[2] = 2.0 * yl;
367: idx[0] = mn - m;
368: idx[1] = mn - m + 1;
369: idx[2] = mn - m - m;
370: PetscCall(MatSetValues(A, 1, &row, 3, idx, v, INSERT_VALUES));
372: row = mn - 1;
373: v[0] = xc;
374: v[1] = 2.0 * xl;
375: v[2] = 2.0 * yl;
376: idx[0] = mn - 1;
377: idx[1] = mn - 2;
378: idx[2] = mn - 1 - m;
379: PetscCall(MatSetValues(A, 1, &i, 3, idx, v, INSERT_VALUES));
381: for (i = 1; i < m - 1; i++) {
382: row = mn - m + i;
383: v[0] = xl;
384: v[1] = xc;
385: v[2] = xr;
386: v[3] = 2.0 * yl;
387: idx[0] = mn - m + i - 1;
388: idx[1] = mn - m + i;
389: idx[2] = mn - m + i + 1;
390: idx[3] = mn - m + i - m;
391: PetscCall(MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES));
392: }
394: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
395: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
396: PetscFunctionReturn(PETSC_SUCCESS);
397: }
399: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
400: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, PetscCtx ctx)
401: {
402: Data *data = (Data *)ctx;
403: PetscInt m, n, mn;
404: PetscReal dx, dy;
405: PetscReal xc, xl, xr, yl, yr;
406: PetscReal a, epsilon;
407: PetscScalar *outptr;
408: const PetscScalar *inptr;
409: PetscInt i, j, len;
410: IS from, to;
411: PetscInt *idx;
412: VecScatter scatter;
413: Vec tmp_in, tmp_out;
415: PetscFunctionBeginUser;
416: m = data->m;
417: n = data->n;
418: mn = m * n;
419: dx = data->dx;
420: dy = data->dy;
421: a = data->a;
422: epsilon = data->epsilon;
424: xc = -2.0 * epsilon * (1.0 / (dx * dx) + 1.0 / (dy * dy));
425: xl = 0.5 * a / dx + epsilon / (dx * dx);
426: xr = -0.5 * a / dx + epsilon / (dx * dx);
427: yl = 0.5 * a / dy + epsilon / (dy * dy);
428: yr = -0.5 * a / dy + epsilon / (dy * dy);
430: /* Get the length of parallel vector */
431: PetscCall(VecGetSize(globalin, &len));
433: /* Set the index sets */
434: PetscCall(PetscMalloc1(len, &idx));
435: for (i = 0; i < len; i++) idx[i] = i;
437: /* Create local sequential vectors */
438: PetscCall(VecCreateSeq(PETSC_COMM_SELF, len, &tmp_in));
439: PetscCall(VecDuplicate(tmp_in, &tmp_out));
441: /* Create scatter context */
442: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, len, idx, PETSC_COPY_VALUES, &from));
443: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, len, idx, PETSC_COPY_VALUES, &to));
444: PetscCall(VecScatterCreate(globalin, from, tmp_in, to, &scatter));
445: PetscCall(VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD));
446: PetscCall(VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD));
447: PetscCall(VecScatterDestroy(&scatter));
449: /*Extract income array - include ghost points */
450: PetscCall(VecGetArrayRead(tmp_in, &inptr));
452: /* Extract outcome array*/
453: PetscCall(VecGetArrayWrite(tmp_out, &outptr));
455: outptr[0] = xc * inptr[0] + xr * inptr[1] + yr * inptr[m];
456: outptr[m - 1] = 2.0 * xl * inptr[m - 2] + xc * inptr[m - 1] + yr * inptr[m - 1 + m];
457: for (i = 1; i < m - 1; i++) outptr[i] = xc * inptr[i] + xl * inptr[i - 1] + xr * inptr[i + 1] + yr * inptr[i + m];
459: for (j = 1; j < n - 1; j++) {
460: outptr[j * m] = xc * inptr[j * m] + xr * inptr[j * m + 1] + yl * inptr[j * m - m] + yr * inptr[j * m + m];
461: 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];
462: 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];
463: }
465: outptr[mn - m] = xc * inptr[mn - m] + xr * inptr[mn - m + 1] + 2.0 * yl * inptr[mn - m - m];
466: outptr[mn - 1] = 2.0 * xl * inptr[mn - 2] + xc * inptr[mn - 1] + 2.0 * yl * inptr[mn - 1 - m];
467: 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];
469: PetscCall(VecRestoreArrayRead(tmp_in, &inptr));
470: PetscCall(VecRestoreArrayWrite(tmp_out, &outptr));
472: PetscCall(VecScatterCreate(tmp_out, from, globalout, to, &scatter));
473: PetscCall(VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));
474: PetscCall(VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));
476: /* Destroy idx and scatter */
477: PetscCall(VecDestroy(&tmp_in));
478: PetscCall(VecDestroy(&tmp_out));
479: PetscCall(ISDestroy(&from));
480: PetscCall(ISDestroy(&to));
481: PetscCall(VecScatterDestroy(&scatter));
483: PetscCall(PetscFree(idx));
484: PetscFunctionReturn(PETSC_SUCCESS);
485: }
487: PetscErrorCode PostStep(TS ts)
488: {
489: PetscReal t;
491: PetscFunctionBeginUser;
492: PetscCall(TSGetTime(ts, &t));
493: PetscCall(PetscPrintf(PETSC_COMM_SELF, " PostStep, t: %g\n", (double)t));
494: PetscFunctionReturn(PETSC_SUCCESS);
495: }
497: /*TEST
499: test:
500: args: -ts_fd -ts_type beuler
501: output_file: output/ex4.out
503: test:
504: suffix: 2
505: args: -ts_fd -ts_type beuler
506: nsize: 2
507: output_file: output/ex4.out
509: test:
510: suffix: 3
511: args: -ts_fd -ts_type cn
513: test:
514: suffix: 4
515: args: -ts_fd -ts_type cn
516: output_file: output/ex4_3.out
517: nsize: 2
519: test:
520: suffix: 5
521: args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
522: output_file: output/ex4.out
524: test:
525: suffix: 6
526: args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
527: output_file: output/ex4.out
528: nsize: 2
530: test:
531: suffix: 7
532: requires: !single
533: args: -ts_fd -ts_type beuler -test_PostStep -ts_time_step .1
535: test:
536: suffix: 8
537: requires: !single
538: args: -ts_type rk -ts_rk_type 5dp -ts_time_step .01 -ts_adapt_type none -ts_view
540: TEST*/