Actual source code: ex4.c
1: /*
2: The Problem:
3: Solve the convection-diffusion equation:
4:
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
10:
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: } Data;
39: int main(int argc,char **argv)
40: {
42: PetscInt time_steps = 100,steps;
43: PetscMPIInt size;
44: Vec global;
45: PetscReal dt,ftime;
46: TS ts;
47: PetscViewer viewfile;
48: MatStructure J_structure;
49: Mat J = 0;
50: Vec x;
51: Data data;
52: PetscInt mn;
53: PetscTruth flg;
54: ISColoring iscoloring;
55: MatFDColoring matfdcoloring = 0;
56: PetscTruth fd_jacobian_coloring = PETSC_FALSE;
57: #if defined(PETSC_HAVE_SUNDIALS)
58: PC pc;
59: PetscViewer viewer;
60: char pcinfo[120],tsinfo[120];
61: const TSType tstype;
62: PetscTruth sundials;
63: #endif
65: PetscInitialize(&argc,&argv,(char*)0,help);
66: MPI_Comm_size(PETSC_COMM_WORLD,&size);
67:
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);
77:
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: TSSetProblemType(ts,TS_NONLINEAR); /* Need to be TS_NONLINEAR for Sundials */
88: TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);
90: /* set user provided RHSFunction and RHSJacobian */
91: TSSetRHSFunction(ts,RHSFunction,&data);
92: MatCreate(PETSC_COMM_WORLD,&J);
93: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
94: MatSetFromOptions(J);
95: MatSeqAIJSetPreallocation(J,5,PETSC_NULL);
96: MatMPIAIJSetPreallocation(J,5,PETSC_NULL,5,PETSC_NULL);
98: PetscOptionsHasName(PETSC_NULL,"-ts_fd",&flg);
99: if (!flg){
100: /* RHSJacobian(ts,0.0,global,&J,&J,&J_structure,&data); */
101: TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
102: } else {
103: PetscOptionsHasName(PETSC_NULL,"-fd_color",&fd_jacobian_coloring);
104: if (fd_jacobian_coloring){ /* Use finite differences with coloring */
105: /* Get data structure of J */
106: PetscTruth pc_diagonal;
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: }
120:
121: /* create coloring context */
122: MatGetColoring(J,MATCOLORING_SL,&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: PetscTruth view_J;
130: PetscOptionsHasName(PETSC_NULL,"-view_J",&view_J);
131: if (view_J){
132: PetscPrintf(PETSC_COMM_SELF,"J computed from TSDefaultComputeJacobianColor():\n");
133: TSDefaultComputeJacobianColor(ts,0.0,global,&J,&J,&J_structure,matfdcoloring);
134: MatView(J,PETSC_VIEWER_STDOUT_WORLD);
135: }
136: } else { /* Use finite differences (slow) */
137: TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobian,&data);
138: }
139: }
141: /* Use SUNDIALS */
142: #if defined(PETSC_HAVE_SUNDIALS)
143: TSSetType(ts,TS_SUNDIALS);
144: #else
145: TSSetType(ts,TS_EULER);
146: #endif
147: dt = 0.1;
148: TSSetInitialTimeStep(ts,0.0,dt);
149: TSSetDuration(ts,time_steps,1);
150: TSSetSolution(ts,global);
151:
152: /* Test TSSetPostStep() and TSSetPostUpdate() */
153: PetscOptionsHasName(PETSC_NULL,"-test_PostStep",&flg);
154: if (flg){
155: TSSetPostStep(ts,PostStep);
156: }
157: PetscOptionsHasName(PETSC_NULL,"-test_PostUpdate",&flg);
158: if (flg){
159: TSSetPostUpdate(ts,PostUpdate);
160: }
162: /* Pick up a Petsc preconditioner */
163: /* one can always set method or preconditioner during the run time */
164: #if defined(PETSC_HAVE_SUNDIALS)
165: TSSundialsGetPC(ts,&pc);
166: PCSetType(pc,PCJACOBI);
167: #endif
168: TSSetFromOptions(ts);
170: TSSetUp(ts);
171: TSStep(ts,&steps,&ftime);
173: PetscOptionsHasName(PETSC_NULL,"-matlab_view",&flg);
174: if (flg){ /* print solution into a MATLAB file */
175: TSGetSolution(ts,&global);
176: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
177: PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
178: VecView(global,viewfile);
179: PetscViewerDestroy(viewfile);
180: }
182: #if defined(PETSC_HAVE_SUNDIALS)
183: /* extracts the PC from ts */
184: TSSundialsGetPC(ts,&pc);
185: TSGetType(ts,&tstype);
186: PetscTypeCompare((PetscObject)ts,TS_SUNDIALS,&sundials);
187: if (sundials){
188: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
189: TSView(ts,viewer);
190: PetscViewerDestroy(viewer);
191: PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
192: PCView(pc,viewer);
193: PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",
194: size,tsinfo,pcinfo);
195: PetscViewerDestroy(viewer);
196: }
197: #endif
199: /* free the memories */
200: TSDestroy(ts);
201: VecDestroy(global);
202: VecDestroy(x);
203: ierr= MatDestroy(J);
204: if (fd_jacobian_coloring){MatFDColoringDestroy(matfdcoloring);}
205: PetscFinalize();
206: return 0;
207: }
209: /* -------------------------------------------------------------------*/
210: /* the initial function */
211: PetscReal f_ini(PetscReal x,PetscReal y)
212: {
213: PetscReal f;
215: f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
216: return f;
217: }
221: PetscErrorCode Initial(Vec global,void *ctx)
222: {
223: Data *data = (Data*)ctx;
224: PetscInt m,row,col;
225: PetscReal x,y,dx,dy;
226: PetscScalar *localptr;
227: PetscInt i,mybase,myend,locsize;
231: /* make the local copies of parameters */
232: m = data->m;
233: dx = data->dx;
234: dy = data->dy;
236: /* determine starting point of each processor */
237: VecGetOwnershipRange(global,&mybase,&myend);
238: VecGetLocalSize(global,&locsize);
240: /* Initialize the array */
241: VecGetArray(global,&localptr);
243: for (i=0; i<locsize; i++) {
244: row = 1+(mybase+i)-((mybase+i)/m)*m;
245: col = (mybase+i)/m+1;
246: x = dx*row;
247: y = dy*col;
248: localptr[i] = f_ini(x,y);
249: }
250:
251: VecRestoreArray(global,&localptr);
252: return(0);
253: }
257: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
258: {
259: VecScatter scatter;
260: IS from,to;
261: PetscInt i,n,*idx,nsteps;
262: Vec tmp_vec;
264: PetscScalar *tmp;
265:
267: TSGetTimeStepNumber(ts,&nsteps);
269: /* Get the size of the vector */
270: VecGetSize(global,&n);
272: /* Set the index sets */
273: PetscMalloc(n*sizeof(PetscInt),&idx);
274: for(i=0; i<n; i++) idx[i]=i;
275:
276: /* Create local sequential vectors */
277: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
279: /* Create scatter context */
280: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
281: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
282: VecScatterCreate(global,from,tmp_vec,to,&scatter);
283: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
284: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
286: VecGetArray(tmp_vec,&tmp);
287: PetscPrintf(PETSC_COMM_WORLD,"At t[%d] =%14.2e u= %14.2e at the center \n",nsteps,time,PetscRealPart(tmp[n/2]));
288: VecRestoreArray(tmp_vec,&tmp);
290: PetscFree(idx);
291: ISDestroy(from);
292: ISDestroy(to);
293: VecScatterDestroy(scatter);
294: VecDestroy(tmp_vec);
295: return(0);
296: }
300: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
301: {
302: Data *data = (Data*)ptr;
303: Mat A = *AA;
304: PetscScalar v[5];
305: PetscInt idx[5],i,j,row;
307: PetscInt m,n,mn;
308: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
311: m = data->m;
312: n = data->n;
313: mn = m*n;
314: dx = data->dx;
315: dy = data->dy;
316: a = data->a;
317: epsilon = data->epsilon;
319: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
320: xl = 0.5*a/dx+epsilon/(dx*dx);
321: xr = -0.5*a/dx+epsilon/(dx*dx);
322: yl = 0.5*a/dy+epsilon/(dy*dy);
323: yr = -0.5*a/dy+epsilon/(dy*dy);
325: row=0;
326: v[0] = xc; v[1]=xr; v[2]=yr;
327: idx[0]=0; idx[1]=2; idx[2]=m;
328: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
330: row=m-1;
331: v[0]=2.0*xl; v[1]=xc; v[2]=yr;
332: idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
333: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
335: for (i=1; i<m-1; i++) {
336: row=i;
337: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
338: idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
339: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
340: }
342: for (j=1; j<n-1; j++) {
343: row=j*m;
344: v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
345: idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
346: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
347:
348: row=j*m+m-1;
349: v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
350: 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;
351: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
353: for(i=1; i<m-1; i++) {
354: row=j*m+i;
355: v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
356: idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
357: idx[4]=j*m+i+m;
358: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
359: }
360: }
362: row=mn-m;
363: v[0] = xc; v[1]=xr; v[2]=2.0*yl;
364: idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
365: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
366:
367: row=mn-1;
368: v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
369: idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
370: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
372: for (i=1; i<m-1; i++) {
373: row=mn-m+i;
374: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
375: idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
376: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
377: }
379: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
380: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
382: /* *flag = SAME_NONZERO_PATTERN; */
383: *flag = DIFFERENT_NONZERO_PATTERN;
384: return(0);
385: }
387: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
390: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
391: {
392: Data *data = (Data*)ctx;
393: PetscInt m,n,mn;
394: PetscReal dx,dy;
395: PetscReal xc,xl,xr,yl,yr;
396: PetscReal a,epsilon;
397: PetscScalar *inptr,*outptr;
398: PetscInt i,j,len;
400: IS from,to;
401: PetscInt *idx;
402: VecScatter scatter;
403: Vec tmp_in,tmp_out;
406: m = data->m;
407: n = data->n;
408: mn = m*n;
409: dx = data->dx;
410: dy = data->dy;
411: a = data->a;
412: epsilon = data->epsilon;
414: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
415: xl = 0.5*a/dx+epsilon/(dx*dx);
416: xr = -0.5*a/dx+epsilon/(dx*dx);
417: yl = 0.5*a/dy+epsilon/(dy*dy);
418: yr = -0.5*a/dy+epsilon/(dy*dy);
419:
420: /* Get the length of parallel vector */
421: VecGetSize(globalin,&len);
423: /* Set the index sets */
424: PetscMalloc(len*sizeof(PetscInt),&idx);
425: for(i=0; i<len; i++) idx[i]=i;
426:
427: /* Create local sequential vectors */
428: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
429: VecDuplicate(tmp_in,&tmp_out);
431: /* Create scatter context */
432: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&from);
433: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&to);
434: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
435: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
436: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
437: VecScatterDestroy(scatter);
439: /*Extract income array - include ghost points */
440: VecGetArray(tmp_in,&inptr);
442: /* Extract outcome array*/
443: VecGetArray(tmp_out,&outptr);
445: outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
446: outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
447: for (i=1; i<m-1; i++) {
448: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
449: +yr*inptr[i+m];
450: }
452: for (j=1; j<n-1; j++) {
453: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
454: yl*inptr[j*m-m]+yr*inptr[j*m+m];
455: outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
456: yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
457: for(i=1; i<m-1; i++) {
458: outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
459: +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
460: }
461: }
463: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
464: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
465: for (i=1; i<m-1; i++) {
466: outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
467: +2*yl*inptr[mn-m+i-m];
468: }
470: VecRestoreArray(tmp_in,&inptr);
471: VecRestoreArray(tmp_out,&outptr);
473: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
474: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
475: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
476:
477: /* Destroy idx aand scatter */
478: VecDestroy(tmp_in);
479: VecDestroy(tmp_out);
480: ISDestroy(from);
481: ISDestroy(to);
482: VecScatterDestroy(scatter);
484: PetscFree(idx);
485: return(0);
486: }
490: PetscErrorCode PostStep(TS ts)
491: {
492: PetscErrorCode ierr;
493: PetscReal t;
494:
496: TSGetTime(ts,&t);
497: PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",t);
498: return(0);
499: }
503: PetscErrorCode PostUpdate(TS ts,PetscReal t, PetscReal *dt)
504: {
505: PetscErrorCode ierr;
507: PetscPrintf(PETSC_COMM_SELF," PostUpdate, t: %g\n",t);
508: return(0);
509: }