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) {
99: PetscCall(TSSetRHSJacobian(ts, J, J, RHSJacobian, &data));
100: } else {
101: PetscCall(TSGetSNES(ts, &snes));
102: PetscCall(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: PetscCall(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: PetscCall(MatGetOwnershipRange(J, &rstart, &rend));
111: for (i = rstart; i < rend; i++) PetscCall(MatSetValues(J, 1, &i, 1, &i, &zero, INSERT_VALUES));
112: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
113: PetscCall(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: PetscCall(TSSetType(ts, TSBEULER));
117: PetscCall(TSSetUp(ts));
118: PetscCall(SNESComputeJacobianDefault(snes, x, J, J, ts));
119: }
121: /* create coloring context */
122: PetscCall(MatColoringCreate(J, &mc));
123: PetscCall(MatColoringSetType(mc, MATCOLORINGSL));
124: PetscCall(MatColoringSetFromOptions(mc));
125: PetscCall(MatColoringApply(mc, &iscoloring));
126: PetscCall(MatColoringDestroy(&mc));
127: PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring));
128: PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode (*)(void))SNESTSFormFunction, ts));
129: PetscCall(MatFDColoringSetFromOptions(matfdcoloring));
130: PetscCall(MatFDColoringSetUp(J, iscoloring, matfdcoloring));
131: PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring));
132: PetscCall(ISColoringDestroy(&iscoloring));
133: } else { /* Use finite differences (slow) */
134: PetscCall(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: PetscCall(TSGetSNES(ts, &snes));
141: PetscCall(SNESGetKSP(snes, &ksp));
142: PetscCall(KSPGetPC(ksp, &pc));
143: PetscCall(PCSetType(pc, PCJACOBI));
144: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
146: PetscCall(TSSetFromOptions(ts));
147: PetscCall(TSSetUp(ts));
149: /* Test TSSetPostStep() */
150: PetscCall(PetscOptionsHasName(NULL, NULL, "-test_PostStep", &flg));
151: if (flg) PetscCall(TSSetPostStep(ts, PostStep));
153: PetscCall(PetscOptionsGetInt(NULL, NULL, "-NOUT", &NOUT, NULL));
154: for (iout = 1; iout <= NOUT; iout++) {
155: PetscCall(TSSetMaxSteps(ts, time_steps));
156: PetscCall(TSSetMaxTime(ts, iout * ftime_original / NOUT));
157: PetscCall(TSSolve(ts, global));
158: PetscCall(TSGetSolveTime(ts, &ftime));
159: PetscCall(TSSetTime(ts, ftime));
160: PetscCall(TSSetTimeStep(ts, dt));
161: }
162: /* Interpolate solution at tfinal */
163: PetscCall(TSGetSolution(ts, &global));
164: PetscCall(TSInterpolate(ts, ftime_original, global));
166: PetscCall(PetscOptionsHasName(NULL, NULL, "-matlab_view", &flg));
167: if (flg) { /* print solution into a MATLAB file */
168: PetscCall(PetscViewerASCIIOpen(PETSC_COMM_WORLD, "out.m", &viewfile));
169: PetscCall(PetscViewerPushFormat(viewfile, PETSC_VIEWER_ASCII_MATLAB));
170: PetscCall(VecView(global, viewfile));
171: PetscCall(PetscViewerPopFormat(viewfile));
172: PetscCall(PetscViewerDestroy(&viewfile));
173: }
175: /* free the memories */
176: PetscCall(TSDestroy(&ts));
177: PetscCall(VecDestroy(&global));
178: PetscCall(VecDestroy(&x));
179: PetscCall(MatDestroy(&J));
180: if (fd_jacobian_coloring) PetscCall(MatFDColoringDestroy(&matfdcoloring));
181: PetscCall(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;
203: PetscFunctionBeginUser;
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: PetscCall(VecGetOwnershipRange(global, &mybase, &myend));
211: PetscCall(VecGetLocalSize(global, &locsize));
213: /* Initialize the array */
214: PetscCall(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: PetscCall(VecRestoreArrayWrite(global, &localptr));
225: PetscFunctionReturn(PETSC_SUCCESS);
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;
236: PetscFunctionBeginUser;
237: PetscCall(TSGetStepNumber(ts, &nsteps));
238: /* display output at selected time steps */
239: PetscCall(TSGetMaxSteps(ts, &maxsteps));
240: if (nsteps % 10 != 0) PetscFunctionReturn(PETSC_SUCCESS);
242: /* Get the size of the vector */
243: PetscCall(VecGetSize(global, &n));
245: /* Set the index sets */
246: PetscCall(PetscMalloc1(n, &idx));
247: for (i = 0; i < n; i++) idx[i] = i;
249: /* Create local sequential vectors */
250: PetscCall(VecCreateSeq(PETSC_COMM_SELF, n, &tmp_vec));
252: /* Create scatter context */
253: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &from));
254: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, n, idx, PETSC_COPY_VALUES, &to));
255: PetscCall(VecScatterCreate(global, from, tmp_vec, to, &scatter));
256: PetscCall(VecScatterBegin(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD));
257: PetscCall(VecScatterEnd(scatter, global, tmp_vec, INSERT_VALUES, SCATTER_FORWARD));
259: PetscCall(VecGetArrayRead(tmp_vec, &tmp));
260: 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])));
261: PetscCall(VecRestoreArrayRead(tmp_vec, &tmp));
263: PetscCall(PetscFree(idx));
264: PetscCall(ISDestroy(&from));
265: PetscCall(ISDestroy(&to));
266: PetscCall(VecScatterDestroy(&scatter));
267: PetscCall(VecDestroy(&tmp_vec));
268: PetscFunctionReturn(PETSC_SUCCESS);
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;
279: PetscFunctionBeginUser;
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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(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: PetscCall(MatSetValues(A, 1, &row, 4, idx, v, INSERT_VALUES));
393: }
395: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
396: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
397: PetscFunctionReturn(PETSC_SUCCESS);
398: }
400: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
401: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec globalin, Vec globalout, void *ctx)
402: {
403: Data *data = (Data *)ctx;
404: PetscInt m, n, mn;
405: PetscReal dx, dy;
406: PetscReal xc, xl, xr, yl, yr;
407: PetscReal a, epsilon;
408: PetscScalar *outptr;
409: const PetscScalar *inptr;
410: PetscInt i, j, len;
411: IS from, to;
412: PetscInt *idx;
413: VecScatter scatter;
414: Vec tmp_in, tmp_out;
416: PetscFunctionBeginUser;
417: m = data->m;
418: n = data->n;
419: mn = m * n;
420: dx = data->dx;
421: dy = data->dy;
422: a = data->a;
423: epsilon = data->epsilon;
425: xc = -2.0 * epsilon * (1.0 / (dx * dx) + 1.0 / (dy * dy));
426: xl = 0.5 * a / dx + epsilon / (dx * dx);
427: xr = -0.5 * a / dx + epsilon / (dx * dx);
428: yl = 0.5 * a / dy + epsilon / (dy * dy);
429: yr = -0.5 * a / dy + epsilon / (dy * dy);
431: /* Get the length of parallel vector */
432: PetscCall(VecGetSize(globalin, &len));
434: /* Set the index sets */
435: PetscCall(PetscMalloc1(len, &idx));
436: for (i = 0; i < len; i++) idx[i] = i;
438: /* Create local sequential vectors */
439: PetscCall(VecCreateSeq(PETSC_COMM_SELF, len, &tmp_in));
440: PetscCall(VecDuplicate(tmp_in, &tmp_out));
442: /* Create scatter context */
443: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, len, idx, PETSC_COPY_VALUES, &from));
444: PetscCall(ISCreateGeneral(PETSC_COMM_SELF, len, idx, PETSC_COPY_VALUES, &to));
445: PetscCall(VecScatterCreate(globalin, from, tmp_in, to, &scatter));
446: PetscCall(VecScatterBegin(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD));
447: PetscCall(VecScatterEnd(scatter, globalin, tmp_in, INSERT_VALUES, SCATTER_FORWARD));
448: PetscCall(VecScatterDestroy(&scatter));
450: /*Extract income array - include ghost points */
451: PetscCall(VecGetArrayRead(tmp_in, &inptr));
453: /* Extract outcome array*/
454: PetscCall(VecGetArrayWrite(tmp_out, &outptr));
456: outptr[0] = xc * inptr[0] + xr * inptr[1] + yr * inptr[m];
457: outptr[m - 1] = 2.0 * xl * inptr[m - 2] + xc * inptr[m - 1] + yr * inptr[m - 1 + m];
458: for (i = 1; i < m - 1; i++) outptr[i] = xc * inptr[i] + xl * inptr[i - 1] + xr * inptr[i + 1] + yr * inptr[i + m];
460: for (j = 1; j < n - 1; j++) {
461: outptr[j * m] = xc * inptr[j * m] + xr * inptr[j * m + 1] + yl * inptr[j * m - m] + yr * inptr[j * m + m];
462: 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];
463: 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];
464: }
466: outptr[mn - m] = xc * inptr[mn - m] + xr * inptr[mn - m + 1] + 2.0 * yl * inptr[mn - m - m];
467: outptr[mn - 1] = 2.0 * xl * inptr[mn - 2] + xc * inptr[mn - 1] + 2.0 * yl * inptr[mn - 1 - m];
468: 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];
470: PetscCall(VecRestoreArrayRead(tmp_in, &inptr));
471: PetscCall(VecRestoreArrayWrite(tmp_out, &outptr));
473: PetscCall(VecScatterCreate(tmp_out, from, globalout, to, &scatter));
474: PetscCall(VecScatterBegin(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));
475: PetscCall(VecScatterEnd(scatter, tmp_out, globalout, INSERT_VALUES, SCATTER_FORWARD));
477: /* Destroy idx and scatter */
478: PetscCall(VecDestroy(&tmp_in));
479: PetscCall(VecDestroy(&tmp_out));
480: PetscCall(ISDestroy(&from));
481: PetscCall(ISDestroy(&to));
482: PetscCall(VecScatterDestroy(&scatter));
484: PetscCall(PetscFree(idx));
485: PetscFunctionReturn(PETSC_SUCCESS);
486: }
488: PetscErrorCode PostStep(TS ts)
489: {
490: PetscReal t;
492: PetscFunctionBeginUser;
493: PetscCall(TSGetTime(ts, &t));
494: PetscCall(PetscPrintf(PETSC_COMM_SELF, " PostStep, t: %g\n", (double)t));
495: PetscFunctionReturn(PETSC_SUCCESS);
496: }
498: /*TEST
500: test:
501: args: -ts_fd -ts_type beuler
502: output_file: output/ex4.out
504: test:
505: suffix: 2
506: args: -ts_fd -ts_type beuler
507: nsize: 2
508: output_file: output/ex4.out
510: test:
511: suffix: 3
512: args: -ts_fd -ts_type cn
514: test:
515: suffix: 4
516: args: -ts_fd -ts_type cn
517: output_file: output/ex4_3.out
518: nsize: 2
520: test:
521: suffix: 5
522: args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
523: output_file: output/ex4.out
525: test:
526: suffix: 6
527: args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
528: output_file: output/ex4.out
529: nsize: 2
531: test:
532: suffix: 7
533: requires: !single
534: args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1
536: test:
537: suffix: 8
538: requires: !single
539: args: -ts_type rk -ts_rk_type 5dp -ts_dt .01 -ts_adapt_type none -ts_view
541: TEST*/