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: {
22: PetscInt m; /* the number of mesh points in x-direction */
23: PetscInt n; /* the number of mesh points in y-direction */
24: PetscReal dx; /* the grid space in x-direction */
25: PetscReal dy; /* the grid space in y-direction */
26: PetscReal a; /* the convection coefficient */
27: PetscReal epsilon; /* the diffusion coefficient */
28: PetscReal tfinal;
29: } Data;
31: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*);
32: extern PetscErrorCode Initial(Vec,void*);
33: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
34: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat,Mat,void*);
35: extern PetscErrorCode PostStep(TS);
37: int main(int argc,char **argv)
38: {
39: PetscInt time_steps=100,iout,NOUT=1;
40: Vec global;
41: PetscReal dt,ftime,ftime_original;
42: TS ts;
43: PetscViewer viewfile;
44: Mat J = 0;
45: Vec x;
46: Data data;
47: PetscInt mn;
48: PetscBool flg;
49: MatColoring mc;
50: ISColoring iscoloring;
51: MatFDColoring matfdcoloring = 0;
52: PetscBool fd_jacobian_coloring = PETSC_FALSE;
53: SNES snes;
54: KSP ksp;
55: 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++) {
112: MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
113: }
114: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
115: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
116: } else {
117: /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
118: TSSetType(ts,TSBEULER);
119: TSSetUp(ts);
120: SNESComputeJacobianDefault(snes,x,J,J,ts);
121: }
123: /* create coloring context */
124: MatColoringCreate(J,&mc);
125: MatColoringSetType(mc,MATCOLORINGSL);
126: MatColoringSetFromOptions(mc);
127: MatColoringApply(mc,&iscoloring);
128: MatColoringDestroy(&mc);
129: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
130: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
131: MatFDColoringSetFromOptions(matfdcoloring);
132: MatFDColoringSetUp(J,iscoloring,matfdcoloring);
133: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
134: ISColoringDestroy(&iscoloring);
135: } else { /* Use finite differences (slow) */
136: SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);
137: }
138: }
140: /* Pick up a Petsc preconditioner */
141: /* one can always set method or preconditioner during the run time */
142: TSGetSNES(ts,&snes);
143: SNESGetKSP(snes,&ksp);
144: KSPGetPC(ksp,&pc);
145: PCSetType(pc,PCJACOBI);
146: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
148: TSSetFromOptions(ts);
149: TSSetUp(ts);
151: /* Test TSSetPostStep() */
152: PetscOptionsHasName(NULL,NULL,"-test_PostStep",&flg);
153: if (flg) {
154: TSSetPostStep(ts,PostStep);
155: }
157: PetscOptionsGetInt(NULL,NULL,"-NOUT",&NOUT,NULL);
158: for (iout=1; iout<=NOUT; iout++) {
159: TSSetMaxSteps(ts,time_steps);
160: TSSetMaxTime(ts,iout*ftime_original/NOUT);
161: TSSolve(ts,global);
162: TSGetSolveTime(ts,&ftime);
163: TSSetTime(ts,ftime);
164: TSSetTimeStep(ts,dt);
165: }
166: /* Interpolate solution at tfinal */
167: TSGetSolution(ts,&global);
168: TSInterpolate(ts,ftime_original,global);
170: PetscOptionsHasName(NULL,NULL,"-matlab_view",&flg);
171: if (flg) { /* print solution into a MATLAB file */
172: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
173: PetscViewerPushFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
174: VecView(global,viewfile);
175: PetscViewerPopFormat(viewfile);
176: PetscViewerDestroy(&viewfile);
177: }
179: /* free the memories */
180: TSDestroy(&ts);
181: VecDestroy(&global);
182: VecDestroy(&x);
183: MatDestroy(&J);
184: if (fd_jacobian_coloring) MatFDColoringDestroy(&matfdcoloring);
185: PetscFinalize();
186: return 0;
187: }
189: /* -------------------------------------------------------------------*/
190: /* the initial function */
191: PetscReal f_ini(PetscReal x,PetscReal y)
192: {
193: PetscReal f;
195: f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2)));
196: return f;
197: }
199: PetscErrorCode Initial(Vec global,void *ctx)
200: {
201: Data *data = (Data*)ctx;
202: PetscInt m,row,col;
203: PetscReal x,y,dx,dy;
204: PetscScalar *localptr;
205: PetscInt i,mybase,myend,locsize;
208: /* make the local copies of parameters */
209: m = data->m;
210: dx = data->dx;
211: dy = data->dy;
213: /* determine starting point of each processor */
214: VecGetOwnershipRange(global,&mybase,&myend);
215: VecGetLocalSize(global,&locsize);
217: /* Initialize the array */
218: VecGetArrayWrite(global,&localptr);
220: for (i=0; i<locsize; i++) {
221: row = 1+(mybase+i)-((mybase+i)/m)*m;
222: col = (mybase+i)/m+1;
223: x = dx*row;
224: y = dy*col;
225: localptr[i] = f_ini(x,y);
226: }
228: VecRestoreArrayWrite(global,&localptr);
229: return 0;
230: }
232: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
233: {
234: VecScatter scatter;
235: IS from,to;
236: PetscInt i,n,*idx,nsteps,maxsteps;
237: Vec tmp_vec;
238: const PetscScalar *tmp;
241: TSGetStepNumber(ts,&nsteps);
242: /* display output at selected time steps */
243: TSGetMaxSteps(ts, &maxsteps);
244: if (nsteps % 10 != 0) return 0;
246: /* Get the size of the vector */
247: VecGetSize(global,&n);
249: /* Set the index sets */
250: PetscMalloc1(n,&idx);
251: for (i=0; i<n; i++) idx[i]=i;
253: /* Create local sequential vectors */
254: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
256: /* Create scatter context */
257: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
258: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
259: VecScatterCreate(global,from,tmp_vec,to,&scatter);
260: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
261: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
263: VecGetArrayRead(tmp_vec,&tmp);
264: PetscPrintf(PETSC_COMM_WORLD,"At t[%D] =%14.2e u= %14.2e at the center \n",nsteps,(double)time,(double)PetscRealPart(tmp[n/2]));
265: VecRestoreArrayRead(tmp_vec,&tmp);
267: PetscFree(idx);
268: ISDestroy(&from);
269: ISDestroy(&to);
270: VecScatterDestroy(&scatter);
271: VecDestroy(&tmp_vec);
272: return 0;
273: }
275: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ptr)
276: {
277: Data *data = (Data*)ptr;
278: PetscScalar v[5];
279: PetscInt idx[5],i,j,row;
280: PetscInt m,n,mn;
281: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
284: m = data->m;
285: n = data->n;
286: mn = m*n;
287: dx = data->dx;
288: dy = data->dy;
289: a = data->a;
290: epsilon = data->epsilon;
292: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
293: xl = 0.5*a/dx+epsilon/(dx*dx);
294: xr = -0.5*a/dx+epsilon/(dx*dx);
295: yl = 0.5*a/dy+epsilon/(dy*dy);
296: yr = -0.5*a/dy+epsilon/(dy*dy);
298: row = 0;
299: v[0] = xc; v[1] = xr; v[2] = yr;
300: idx[0] = 0; idx[1] = 2; idx[2] = m;
301: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
303: row = m-1;
304: v[0] = 2.0*xl; v[1] = xc; v[2] = yr;
305: idx[0] = m-2; idx[1] = m-1; idx[2] = m-1+m;
306: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
308: for (i=1; i<m-1; i++) {
309: row = i;
310: v[0] = xl; v[1] = xc; v[2] = xr; v[3] = yr;
311: idx[0] = i-1; idx[1] = i; idx[2] = i+1; idx[3] = i+m;
312: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
313: }
315: for (j=1; j<n-1; j++) {
316: row = j*m;
317: v[0] = xc; v[1] = xr; v[2] = yl; v[3] = yr;
318: idx[0] = j*m; idx[1] = j*m; idx[2] = j*m-m; idx[3] = j*m+m;
319: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
321: row = j*m+m-1;
322: v[0] = xc; v[1] = 2.0*xl; v[2] = yl; v[3] = yr;
323: idx[0] = j*m+m-1; idx[1] = j*m+m-1-1; idx[2] = j*m+m-1-m; idx[3] = j*m+m-1+m;
324: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
326: for (i=1; i<m-1; i++) {
327: row = j*m+i;
328: v[0] = xc; v[1] = xl; v[2] = xr; v[3] = yl; v[4]=yr;
329: idx[0] = j*m+i; idx[1] = j*m+i-1; idx[2] = j*m+i+1; idx[3] = j*m+i-m;
330: idx[4] = j*m+i+m;
331: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
332: }
333: }
335: row = mn-m;
336: v[0] = xc; v[1] = xr; v[2] = 2.0*yl;
337: idx[0] = mn-m; idx[1] = mn-m+1; idx[2] = mn-m-m;
338: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
340: row = mn-1;
341: v[0] = xc; v[1] = 2.0*xl; v[2] = 2.0*yl;
342: idx[0] = mn-1; idx[1] = mn-2; idx[2] = mn-1-m;
343: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
345: for (i=1; i<m-1; i++) {
346: row = mn-m+i;
347: v[0] = xl; v[1] = xc; v[2] = xr; v[3] = 2.0*yl;
348: idx[0] = mn-m+i-1; idx[1] = mn-m+i; idx[2] = mn-m+i+1; idx[3] = mn-m+i-m;
349: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
350: }
352: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
353: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
355: return 0;
356: }
358: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
359: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
360: {
361: Data *data = (Data*)ctx;
362: PetscInt m,n,mn;
363: PetscReal dx,dy;
364: PetscReal xc,xl,xr,yl,yr;
365: PetscReal a,epsilon;
366: PetscScalar *outptr;
367: const PetscScalar *inptr;
368: PetscInt i,j,len;
369: IS from,to;
370: PetscInt *idx;
371: VecScatter scatter;
372: Vec tmp_in,tmp_out;
375: m = data->m;
376: n = data->n;
377: mn = m*n;
378: dx = data->dx;
379: dy = data->dy;
380: a = data->a;
381: epsilon = data->epsilon;
383: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
384: xl = 0.5*a/dx+epsilon/(dx*dx);
385: xr = -0.5*a/dx+epsilon/(dx*dx);
386: yl = 0.5*a/dy+epsilon/(dy*dy);
387: yr = -0.5*a/dy+epsilon/(dy*dy);
389: /* Get the length of parallel vector */
390: VecGetSize(globalin,&len);
392: /* Set the index sets */
393: PetscMalloc1(len,&idx);
394: for (i=0; i<len; i++) idx[i]=i;
396: /* Create local sequential vectors */
397: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
398: VecDuplicate(tmp_in,&tmp_out);
400: /* Create scatter context */
401: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
402: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
403: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
404: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
405: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
406: VecScatterDestroy(&scatter);
408: /*Extract income array - include ghost points */
409: VecGetArrayRead(tmp_in,&inptr);
411: /* Extract outcome array*/
412: VecGetArrayWrite(tmp_out,&outptr);
414: outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
415: outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
416: for (i=1; i<m-1; i++) {
417: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]+yr*inptr[i+m];
418: }
420: for (j=1; j<n-1; j++) {
421: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+ yl*inptr[j*m-m]+yr*inptr[j*m+m];
422: 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];
423: for (i=1; i<m-1; i++) {
424: 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];
425: }
426: }
428: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
429: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
430: for (i=1; i<m-1; i++) {
431: 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];
432: }
434: VecRestoreArrayRead(tmp_in,&inptr);
435: VecRestoreArrayWrite(tmp_out,&outptr);
437: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
438: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
439: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
441: /* Destroy idx aand scatter */
442: VecDestroy(&tmp_in);
443: VecDestroy(&tmp_out);
444: ISDestroy(&from);
445: ISDestroy(&to);
446: VecScatterDestroy(&scatter);
448: PetscFree(idx);
449: return 0;
450: }
452: PetscErrorCode PostStep(TS ts)
453: {
454: PetscReal t;
457: TSGetTime(ts,&t);
458: PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",(double)t);
459: return 0;
460: }
462: /*TEST
464: test:
465: args: -ts_fd -ts_type beuler
466: output_file: output/ex4.out
468: test:
469: suffix: 2
470: args: -ts_fd -ts_type beuler
471: nsize: 2
472: output_file: output/ex4.out
474: test:
475: suffix: 3
476: args: -ts_fd -ts_type cn
478: test:
479: suffix: 4
480: args: -ts_fd -ts_type cn
481: output_file: output/ex4_3.out
482: nsize: 2
484: test:
485: suffix: 5
486: args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
487: output_file: output/ex4.out
489: test:
490: suffix: 6
491: args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
492: output_file: output/ex4.out
493: nsize: 2
495: test:
496: suffix: 7
497: requires: !single
498: args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1
500: test:
501: suffix: 8
502: requires: !single
503: args: -ts_type rk -ts_rk_type 5dp -ts_dt .01 -ts_adapt_type none -ts_view
505: TEST*/