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 with SUNDIALS.
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;
39: int main(int argc,char **argv)
40: {
42: PetscInt time_steps=100,iout,NOUT=1;
43: PetscMPIInt size;
44: Vec global;
45: PetscReal dt,ftime,ftime_original;
46: TS ts;
47: PetscViewer viewfile;
48: MatStructure J_structure;
49: Mat J = 0;
50: Vec x;
51: Data data;
52: PetscInt mn;
53: PetscBool flg;
54: ISColoring iscoloring;
55: MatFDColoring matfdcoloring = 0;
56: PetscBool fd_jacobian_coloring = PETSC_FALSE;
57: SNES snes;
58: KSP ksp;
59: PC pc;
60: PetscViewer viewer;
61: char pcinfo[120],tsinfo[120];
62: const TSType tstype;
63: PetscBool sundials;
65: PetscInitialize(&argc,&argv,(char*)0,help);
66: MPI_Comm_size(PETSC_COMM_WORLD,&size);
68: /* set data */
69: data.m = 9;
70: data.n = 9;
71: data.a = 1.0;
72: data.epsilon = 0.1;
73: data.dx = 1.0/(data.m+1.0);
74: data.dy = 1.0/(data.n+1.0);
75: mn = (data.m)*(data.n);
76: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
78: /* set initial conditions */
79: VecCreate(PETSC_COMM_WORLD,&global);
80: VecSetSizes(global,PETSC_DECIDE,mn);
81: VecSetFromOptions(global);
82: Initial(global,&data);
83: VecDuplicate(global,&x);
85: /* create timestep context */
86: TSCreate(PETSC_COMM_WORLD,&ts);
87: TSMonitorSet(ts,Monitor,&data,PETSC_NULL);
89: /* set user provided RHSFunction and RHSJacobian */
90: TSSetRHSFunction(ts,PETSC_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,PETSC_NULL);
95: MatMPIAIJSetPreallocation(J,5,PETSC_NULL,5,PETSC_NULL);
97: PetscOptionsHasName(PETSC_NULL,"-ts_fd",&flg);
98: if (!flg){
99: /* RHSJacobian(ts,0.0,global,&J,&J,&J_structure,&data); */
100: TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
101: } else {
102: PetscOptionsHasName(PETSC_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: PetscBool view_J;
107: PetscOptionsHasName(PETSC_NULL,"-pc_diagonal",&pc_diagonal);
108: if (pc_diagonal){ /* the preconditioner of J is a diagonal matrix */
109: PetscInt rstart,rend,i;
110: PetscScalar zero=0.0;
111: MatGetOwnershipRange(J,&rstart,&rend);
112: for (i=rstart; i<rend; i++){
113: MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
114: }
115: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
116: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
117: } else { /* the preconditioner of J has same data structure as J */
118: TSDefaultComputeJacobian(ts,0.0,global,&J,&J,&J_structure,&data);
119: }
121: /* create coloring context */
122: MatGetColoring(J,MATCOLORINGSL,&iscoloring);
123: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
124: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))RHSFunction,&data);
125: MatFDColoringSetFromOptions(matfdcoloring);
126: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
127: ISColoringDestroy(&iscoloring);
129: PetscOptionsHasName(PETSC_NULL,"-view_J",&view_J);
130: if (view_J){
131: PetscPrintf(PETSC_COMM_SELF,"J computed from TSDefaultComputeJacobianColor():\n");
132: TSDefaultComputeJacobianColor(ts,0.0,global,&J,&J,&J_structure,matfdcoloring);
133: MatView(J,PETSC_VIEWER_STDOUT_WORLD);
134: }
135: } else { /* Use finite differences (slow) */
136: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobian,&data);
137: }
138: }
140: /* Use SUNDIALS */
141: #if defined(PETSC_HAVE_SUNDIALS)
142: TSSetType(ts,TSSUNDIALS);
143: #else
144: TSSetType(ts,TSEULER);
145: #endif
146: dt = 0.1;
147: ftime_original = data.tfinal = 1.0;
148: TSSetInitialTimeStep(ts,0.0,dt);
149: TSSetDuration(ts,time_steps,ftime_original);
150: TSSetSolution(ts,global);
153: /* Pick up a Petsc preconditioner */
154: /* one can always set method or preconditioner during the run time */
155: TSGetSNES(ts,&snes);
156: SNESGetKSP(snes,&ksp);
157: KSPGetPC(ksp,&pc);
158: PCSetType(pc,PCJACOBI);
160: TSSetFromOptions(ts);
161: TSSetUp(ts);
162:
163: /* Test TSSetPostStep() */
164: PetscOptionsHasName(PETSC_NULL,"-test_PostStep",&flg);
165: if (flg){
166: TSSetPostStep(ts,PostStep);
167: }
169: PetscOptionsGetInt(PETSC_NULL,"-NOUT",&NOUT,PETSC_NULL);
170: for (iout=1; iout<=NOUT; iout++){
171: TSSetDuration(ts,time_steps,iout*ftime_original/NOUT);
172: TSSolve(ts,global,&ftime);
173: TSSetInitialTimeStep(ts,ftime,dt);
174: }
175: /* Interpolate solution at tfinal */
176: TSGetSolution(ts,&global);
177: TSInterpolate(ts,ftime_original,global);
179: PetscOptionsHasName(PETSC_NULL,"-matlab_view",&flg);
180: if (flg){ /* print solution into a MATLAB file */
181: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
182: PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
183: VecView(global,viewfile);
184: PetscViewerDestroy(&viewfile);
185: }
187: /* display solver info for Sundials */
188: TSGetType(ts,&tstype);
189: PetscTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);
190: if (sundials){
191: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
192: TSView(ts,viewer);
193: PetscViewerDestroy(&viewer);
194: PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
195: PCView(pc,viewer);
196: PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",
197: size,tsinfo,pcinfo);
198: PetscViewerDestroy(&viewer);
199: }
201: /* free the memories */
202: TSDestroy(&ts);
203: VecDestroy(&global);
204: VecDestroy(&x);
205: ierr= MatDestroy(&J);
206: if (fd_jacobian_coloring){MatFDColoringDestroy(&matfdcoloring);}
207: PetscFinalize();
208: return 0;
209: }
211: /* -------------------------------------------------------------------*/
212: /* the initial function */
213: PetscReal f_ini(PetscReal x,PetscReal y)
214: {
215: PetscReal f;
217: f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
218: return f;
219: }
223: PetscErrorCode Initial(Vec global,void *ctx)
224: {
225: Data *data = (Data*)ctx;
226: PetscInt m,row,col;
227: PetscReal x,y,dx,dy;
228: PetscScalar *localptr;
229: PetscInt i,mybase,myend,locsize;
233: /* make the local copies of parameters */
234: m = data->m;
235: dx = data->dx;
236: dy = data->dy;
238: /* determine starting point of each processor */
239: VecGetOwnershipRange(global,&mybase,&myend);
240: VecGetLocalSize(global,&locsize);
242: /* Initialize the array */
243: VecGetArray(global,&localptr);
245: for (i=0; i<locsize; i++) {
246: row = 1+(mybase+i)-((mybase+i)/m)*m;
247: col = (mybase+i)/m+1;
248: x = dx*row;
249: y = dy*col;
250: localptr[i] = f_ini(x,y);
251: }
253: VecRestoreArray(global,&localptr);
254: return(0);
255: }
259: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
260: {
261: VecScatter scatter;
262: IS from,to;
263: PetscInt i,n,*idx,nsteps,maxsteps;
264: Vec tmp_vec;
266: PetscScalar *tmp;
267: PetscReal maxtime;
268: Data *data = (Data*)ctx;
269: PetscReal tfinal = data->tfinal;
272: if (time > tfinal) return(0);
274: TSGetTimeStepNumber(ts,&nsteps);
275: /* display output at selected time steps */
276: TSGetDuration(ts, &maxsteps, &maxtime);
277: if (nsteps % 10 != 0 && time < maxtime) return(0);
279: /* Get the size of the vector */
280: VecGetSize(global,&n);
282: /* Set the index sets */
283: PetscMalloc(n*sizeof(PetscInt),&idx);
284: for(i=0; i<n; i++) idx[i]=i;
286: /* Create local sequential vectors */
287: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
289: /* Create scatter context */
290: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
291: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
292: VecScatterCreate(global,from,tmp_vec,to,&scatter);
293: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
294: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
296: VecGetArray(tmp_vec,&tmp);
297: PetscPrintf(PETSC_COMM_WORLD,"At t[%d] =%14.2e u= %14.2e at the center \n",nsteps,time,PetscRealPart(tmp[n/2]));
298: VecRestoreArray(tmp_vec,&tmp);
300: PetscFree(idx);
301: ISDestroy(&from);
302: ISDestroy(&to);
303: VecScatterDestroy(&scatter);
304: VecDestroy(&tmp_vec);
305: return(0);
306: }
310: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
311: {
312: Data *data = (Data*)ptr;
313: Mat A = *AA;
314: PetscScalar v[5];
315: PetscInt idx[5],i,j,row;
317: PetscInt m,n,mn;
318: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
321: m = data->m;
322: n = data->n;
323: mn = m*n;
324: dx = data->dx;
325: dy = data->dy;
326: a = data->a;
327: epsilon = data->epsilon;
329: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
330: xl = 0.5*a/dx+epsilon/(dx*dx);
331: xr = -0.5*a/dx+epsilon/(dx*dx);
332: yl = 0.5*a/dy+epsilon/(dy*dy);
333: yr = -0.5*a/dy+epsilon/(dy*dy);
335: row=0;
336: v[0] = xc; v[1]=xr; v[2]=yr;
337: idx[0]=0; idx[1]=2; idx[2]=m;
338: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
340: row=m-1;
341: v[0]=2.0*xl; v[1]=xc; v[2]=yr;
342: idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
343: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
345: for (i=1; i<m-1; i++) {
346: row=i;
347: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
348: idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
349: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
350: }
352: for (j=1; j<n-1; j++) {
353: row=j*m;
354: v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
355: idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
356: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
358: row=j*m+m-1;
359: v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
360: 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;
361: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
363: for(i=1; i<m-1; i++) {
364: row=j*m+i;
365: v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
366: idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
367: idx[4]=j*m+i+m;
368: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
369: }
370: }
372: row=mn-m;
373: v[0] = xc; v[1]=xr; v[2]=2.0*yl;
374: idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
375: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
377: row=mn-1;
378: v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
379: idx[0]=mn-1; idx[1]=mn-2; 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; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
385: idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
386: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
387: }
389: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
390: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
392: /* *flag = SAME_NONZERO_PATTERN; */
393: *flag = DIFFERENT_NONZERO_PATTERN;
394: return(0);
395: }
397: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
400: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *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 *inptr,*outptr;
408: PetscInt i,j,len;
410: IS from,to;
411: PetscInt *idx;
412: VecScatter scatter;
413: Vec tmp_in,tmp_out;
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: VecGetSize(globalin,&len);
433: /* Set the index sets */
434: PetscMalloc(len*sizeof(PetscInt),&idx);
435: for(i=0; i<len; i++) idx[i]=i;
437: /* Create local sequential vectors */
438: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
439: VecDuplicate(tmp_in,&tmp_out);
441: /* Create scatter context */
442: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
443: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
444: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
445: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
446: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
447: VecScatterDestroy(&scatter);
449: /*Extract income array - include ghost points */
450: VecGetArray(tmp_in,&inptr);
452: /* Extract outcome array*/
453: VecGetArray(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++) {
458: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
459: +yr*inptr[i+m];
460: }
462: for (j=1; j<n-1; j++) {
463: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
464: yl*inptr[j*m-m]+yr*inptr[j*m+m];
465: outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
466: yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
467: for(i=1; i<m-1; i++) {
468: outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
469: +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
470: }
471: }
473: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
474: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
475: for (i=1; i<m-1; i++) {
476: outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
477: +2*yl*inptr[mn-m+i-m];
478: }
480: VecRestoreArray(tmp_in,&inptr);
481: VecRestoreArray(tmp_out,&outptr);
483: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
484: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
485: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
487: /* Destroy idx aand scatter */
488: VecDestroy(&tmp_in);
489: VecDestroy(&tmp_out);
490: ISDestroy(&from);
491: ISDestroy(&to);
492: VecScatterDestroy(&scatter);
494: PetscFree(idx);
495: return(0);
496: }
500: PetscErrorCode PostStep(TS ts)
501: {
502: PetscErrorCode ierr;
503: PetscReal t;
506: TSGetTime(ts,&t);
507: PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",t);
508: return(0);
509: }