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 petscsys.h
19: #include petscts.h
24: typedef struct
25: {
26: PetscInt m; /* the number of mesh points in x-direction */
27: PetscInt n; /* the number of mesh points in y-direction */
28: PetscReal dx; /* the grid space in x-direction */
29: PetscReal dy; /* the grid space in y-direction */
30: PetscReal a; /* the convection coefficient */
31: PetscReal epsilon; /* the diffusion coefficient */
32: } Data;
34: /* two temporal functions */
40: /* #undef PETSC_HAVE_SUNDIALS */
44: int main(int argc,char **argv)
45: {
47: PetscInt time_steps = 100,steps;
48: PetscMPIInt size;
49: Vec global;
50: PetscReal dt,ftime;
51: TS ts;
52: PetscViewer viewfile;
53: MatStructure J_structure;
54: Mat J = 0;
55: //SNES snes;
56: Vec x;
57: Data data;
58: PetscInt mn;
59: PetscTruth flg;
60: #if defined(PETSC_HAVE_SUNDIALS)
61: PC pc;
62: PetscViewer viewer;
63: char pcinfo[120],tsinfo[120];
64: #endif
66: PetscInitialize(&argc,&argv,(char*)0,help);
67: MPI_Comm_size(PETSC_COMM_WORLD,&size);
68:
69: /* set data */
70: data.m = 9;
71: data.n = 9;
72: data.a = 1.0;
73: data.epsilon = 0.1;
74: data.dx = 1.0/(data.m+1.0);
75: data.dy = 1.0/(data.n+1.0);
76: mn = (data.m)*(data.n);
77: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
78:
79: /* set initial conditions */
80: VecCreate(PETSC_COMM_WORLD,&global);
81: VecSetSizes(global,PETSC_DECIDE,mn);
82: VecSetFromOptions(global);
83: Initial(global,&data);
84: VecDuplicate(global,&x);
86: /* create timestep context */
87: TSCreate(PETSC_COMM_WORLD,&ts);
88: TSSetProblemType(ts,TS_NONLINEAR); /* Need to be TS_NONLINEAR for Sundials */
89: TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);
91: /* set user provided RHSFunction and RHSJacobian */
92: TSSetRHSFunction(ts,RHSFunction,&data);
93: MatCreate(PETSC_COMM_WORLD,&J);
94: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
95: MatSetFromOptions(J);
96: RHSJacobian(ts,0.0,global,&J,&J,&J_structure,&data);
97: TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
99: /* Use SUNDIALS */
100: #if defined(PETSC_HAVE_SUNDIALS)
101: TSSetType(ts,TS_SUNDIALS);
102: #else
103: TSSetType(ts,TS_EULER);
104: #endif
105: dt = 0.1;
106: TSSetInitialTimeStep(ts,0.0,dt);
107: TSSetDuration(ts,time_steps,1);
108: TSSetSolution(ts,global);
111: /* Pick up a Petsc preconditioner */
112: /* one can always set method or preconditioner during the run time */
113: #if defined(PETSC_HAVE_SUNDIALS)
114: TSSundialsGetPC(ts,&pc);
115: PCSetType(pc,PCJACOBI);
116: //TSSundialsSetType(ts,SUNDIALS_ADAMS);
117: #endif
118: TSSetFromOptions(ts);
120: TSSetUp(ts);
121: TSStep(ts,&steps,&ftime);
123: PetscOptionsHasName(PETSC_NULL,"-matlab_view",&flg);
124: if (flg){ /* print solution into a MATLAB file */
125: TSGetSolution(ts,&global);
126: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
127: PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
128: VecView(global,viewfile);
129: PetscViewerDestroy(viewfile);
130: }
132: #if defined(PETSC_HAVE_SUNDIALS)
133: /* extracts the PC from ts */
134: TSSundialsGetPC(ts,&pc);
135: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
136: TSView(ts,viewer);
137: PetscViewerDestroy(viewer);
138: PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
139: PCView(pc,viewer);
140: PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",
141: size,tsinfo,pcinfo);
142: PetscViewerDestroy(viewer);
143: #endif
145: /* free the memories */
146: TSDestroy(ts);
147: VecDestroy(global);
148: VecDestroy(x);
149: if (J) {ierr= MatDestroy(J);}
150: //SNESDestroy(snes);
151: PetscFinalize();
152: return 0;
153: }
155: /* -------------------------------------------------------------------*/
156: /* the initial function */
157: PetscReal f_ini(PetscReal x,PetscReal y)
158: {
159: PetscReal f;
160: f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
161: return f;
162: }
166: PetscErrorCode Initial(Vec global,void *ctx)
167: {
168: Data *data = (Data*)ctx;
169: PetscInt m,row,col;
170: PetscReal x,y,dx,dy;
171: PetscScalar *localptr;
172: PetscInt i,mybase,myend,locsize;
175: /* make the local copies of parameters */
176: m = data->m;
177: dx = data->dx;
178: dy = data->dy;
180: /* determine starting point of each processor */
181: VecGetOwnershipRange(global,&mybase,&myend);
182: VecGetLocalSize(global,&locsize);
184: /* Initialize the array */
185: VecGetArray(global,&localptr);
187: for (i=0; i<locsize; i++) {
188: row = 1+(mybase+i)-((mybase+i)/m)*m;
189: col = (mybase+i)/m+1;
190: x = dx*row;
191: y = dy*col;
192: localptr[i] = f_ini(x,y);
193: }
194:
195: VecRestoreArray(global,&localptr);
196: return 0;
197: }
201: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
202: {
203: VecScatter scatter;
204: IS from,to;
205: PetscInt i,n,*idx;
206: Vec tmp_vec;
208: PetscScalar *tmp;
210: /* Get the size of the vector */
211: VecGetSize(global,&n);
213: /* Set the index sets */
214: PetscMalloc(n*sizeof(PetscInt),&idx);
215: for(i=0; i<n; i++) idx[i]=i;
216:
217: /* Create local sequential vectors */
218: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
220: /* Create scatter context */
221: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
222: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
223: VecScatterCreate(global,from,tmp_vec,to,&scatter);
224: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
225: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
227: VecGetArray(tmp_vec,&tmp);
228: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u= %14.6e at the center \n",time,PetscRealPart(tmp[n/2]));
229: VecRestoreArray(tmp_vec,&tmp);
231: PetscFree(idx);
232: ISDestroy(from);
233: ISDestroy(to);
234: VecScatterDestroy(scatter);
235: VecDestroy(tmp_vec);
236: return 0;
237: }
239: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
242: PetscErrorCode FormFunction(SNES snes,Vec globalin,Vec globalout,void *ptr)
243: {
244: Data *data = (Data*)ptr;
245: PetscInt m,n,mn;
246: PetscReal dx,dy;
247: PetscReal xc,xl,xr,yl,yr;
248: PetscReal a,epsilon;
249: PetscScalar *inptr,*outptr;
250: PetscInt i,j,len;
252: IS from,to;
253: PetscInt *idx;
254: VecScatter scatter;
255: Vec tmp_in,tmp_out;
257: m = data->m;
258: n = data->n;
259: mn = m*n;
260: dx = data->dx;
261: dy = data->dy;
262: a = data->a;
263: epsilon = data->epsilon;
265: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
266: xl = 0.5*a/dx+epsilon/(dx*dx);
267: xr = -0.5*a/dx+epsilon/(dx*dx);
268: yl = 0.5*a/dy+epsilon/(dy*dy);
269: yr = -0.5*a/dy+epsilon/(dy*dy);
270:
271: /* Get the length of parallel vector */
272: VecGetSize(globalin,&len);
274: /* Set the index sets */
275: PetscMalloc(len*sizeof(PetscInt),&idx);
276: for(i=0; i<len; i++) idx[i]=i;
277:
278: /* Create local sequential vectors */
279: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
280: VecDuplicate(tmp_in,&tmp_out);
282: /* Create scatter context */
283: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&from);
284: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&to);
285: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
286: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
287: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
288: VecScatterDestroy(scatter);
290: /*Extract income array - include ghost points */
291: VecGetArray(tmp_in,&inptr);
293: /* Extract outcome array*/
294: VecGetArray(tmp_out,&outptr);
296: outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
297: outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
298: for (i=1; i<m-1; i++) {
299: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
300: +yr*inptr[i+m];
301: }
303: for (j=1; j<n-1; j++) {
304: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
305: yl*inptr[j*m-m]+yr*inptr[j*m+m];
306: outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
307: yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
308: for(i=1; i<m-1; i++) {
309: outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
310: +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
311: }
312: }
314: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
315: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
316: for (i=1; i<m-1; i++) {
317: outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
318: +2*yl*inptr[mn-m+i-m];
319: }
321: VecRestoreArray(tmp_in,&inptr);
322: VecRestoreArray(tmp_out,&outptr);
324: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
325: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
326: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
327:
328: /* Destroy idx aand scatter */
329: VecDestroy(tmp_in);
330: VecDestroy(tmp_out);
331: ISDestroy(from);
332: ISDestroy(to);
333: VecScatterDestroy(scatter);
335: PetscFree(idx);
336: return 0;
337: }
341: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
342: {
343: Data *data = (Data*)ptr;
344: Mat A = *AA;
345: PetscScalar v[1],one = 1.0;
346: PetscInt idx[1],i,j,row;
348: PetscInt m,n,mn;
350: m = data->m;
351: n = data->n;
352: mn = (data->m)*(data->n);
353:
354: for(i=0; i<mn; i++) {
355: idx[0] = i;
356: v[0] = one;
357: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
358: }
360: for(i=0; i<mn-m; i++) {
361: idx[0] = i+m;
362: v[0]= one;
363: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
364: }
366: for(i=m; i<mn; i++) {
367: idx[0] = i-m;
368: v[0] = one;
369: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
370: }
372: for(j=0; j<n; j++) {
373: for(i=0; i<m-1; i++) {
374: row = j*m+i;
375: idx[0] = j*m+i+1;
376: v[0] = one;
377: MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
378: }
379: }
381: for(j=0; j<n; j++) {
382: for(i=1; i<m; i++) {
383: row = j*m+i;
384: idx[0] = j*m+i-1;
385: v[0] = one;
386: MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
387: }
388: }
389: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
390: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
391: *flag = SAME_NONZERO_PATTERN;
392: return 0;
393: }
397: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
398: {
399: Data *data = (Data*)ptr;
400: Mat A = *AA;
401: PetscScalar v[5];
402: PetscInt idx[5],i,j,row;
404: PetscInt m,n,mn;
405: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
407: m = data->m;
408: n = data->n;
409: mn = m*n;
410: dx = data->dx;
411: dy = data->dy;
412: a = data->a;
413: epsilon = data->epsilon;
415: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
416: xl = 0.5*a/dx+epsilon/(dx*dx);
417: xr = -0.5*a/dx+epsilon/(dx*dx);
418: yl = 0.5*a/dy+epsilon/(dy*dy);
419: yr = -0.5*a/dy+epsilon/(dy*dy);
421: row=0;
422: v[0] = xc; v[1]=xr; v[2]=yr;
423: idx[0]=0; idx[1]=2; idx[2]=m;
424: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
426: row=m-1;
427: v[0]=2.0*xl; v[1]=xc; v[2]=yr;
428: idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
429: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
431: for (i=1; i<m-1; i++) {
432: row=i;
433: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
434: idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
435: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
436: }
438: for (j=1; j<n-1; j++) {
439: row=j*m;
440: v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
441: idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
442: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
443:
444: row=j*m+m-1;
445: v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
446: 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;
447: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
449: for(i=1; i<m-1; i++) {
450: row=j*m+i;
451: v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
452: idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
453: idx[4]=j*m+i+m;
454: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
455: }
456: }
458: row=mn-m;
459: v[0] = xc; v[1]=xr; v[2]=2.0*yl;
460: idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
461: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
462:
463: row=mn-1;
464: v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
465: idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
466: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
468: for (i=1; i<m-1; i++) {
469: row=mn-m+i;
470: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
471: idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
472: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
473: }
476: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
477: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
479: /* *flag = SAME_NONZERO_PATTERN; */
480: *flag = DIFFERENT_NONZERO_PATTERN;
481: return 0;
482: }
486: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
487: {
489: SNES snes = PETSC_NULL;
491: FormFunction(snes,globalin,globalout,ctx);
492: return 0;
493: }