Actual source code: ex4.c
1: /*$Id: ex4.c,v 1.13 2001/08/07 21:31:27 bsmith Exp $*/
2: /*
3: The Problem:
4: Solve the convection-diffusion equation:
5:
6: u_t+a*(u_x+u_y)=epsilon*(u_xx+u_yy)
7: u=0 at x=0, y=0
8: u_x=0 at x=1
9: u_y=0 at y=1
10:
11: This program tests the routine of computing the Jacobian by the
12: finite difference method as well as PETSc with PVODE.
14: */
16: static char help[] = "Solve the convection-diffusion equation. \n\n";
18: #include petscsys.h
19: #include petscts.h
21: extern int Monitor(TS,int,PetscReal,Vec,void *);
22: extern int Initial(Vec,void *);
24: typedef struct
25: {
26: int m; /* the number of mesh points in x-direction */
27: int 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 */
35: extern int FormJacobian(SNES,Vec,Mat*,Mat*,MatStructure*,void*);
36: extern int FormFunction(SNES,Vec,Vec,void*);
37: extern int RHSFunction(TS,PetscReal,Vec,Vec,void*);
38: extern int RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure *,void*);
40: /* the initial function */
41: PetscReal f_ini(PetscReal x,PetscReal y)
42: {
43: PetscReal f;
44: f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
45: return f;
46: }
49: #define linear_no_matrix 0
50: #define linear_no_time 1
51: #define linear 2
52: #define nonlinear_no_jacobian 3
53: #define nonlinear 4
57: int main(int argc,char **argv)
58: {
59: int ierr,time_steps = 100,steps,size;
60: Vec global;
61: PetscReal dt,ftime;
62: TS ts;
63: PetscViewer viewfile;
64: MatStructure A_structure;
65: Mat A = 0;
66: TSProblemType tsproblem = TS_NONLINEAR; /* Need to be TS_NONLINEAR */
67: SNES snes;
68: Vec x;
69: Data data;
70: int mn;
71: #if defined(PETSC_HAVE_PVODE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
72: PC pc;
73: PetscViewer viewer;
74: char pcinfo[120],tsinfo[120];
75: #endif
77: PetscInitialize(&argc,&argv,(char*)0,help);
78: MPI_Comm_size(PETSC_COMM_WORLD,&size);
79:
80: /* set Data */
81: data.m = 9;
82: data.n = 9;
83: data.a = 1.0;
84: data.epsilon = 0.1;
85: data.dx = 1.0/(data.m+1.0);
86: data.dy = 1.0/(data.n+1.0);
87: mn = (data.m)*(data.n);
89: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
90:
91: /* set initial conditions */
92: VecCreate(PETSC_COMM_WORLD,&global);
93: VecSetSizes(global,PETSC_DECIDE,mn);
94: VecSetFromOptions(global);
95: Initial(global,&data);
96: VecDuplicate(global,&x);
97:
98: /* make timestep context */
99: TSCreate(PETSC_COMM_WORLD,&ts);
100: TSSetProblemType(ts,tsproblem);
101: TSSetMonitor(ts,Monitor,PETSC_NULL,PETSC_NULL);
103: dt = 0.1;
105: /*
106: The user provides the RHS and Jacobian
107: */
108: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,mn,mn,&A);
109: MatSetFromOptions(A);
111: /* Create SNES context */
112: SNESCreate(PETSC_COMM_WORLD,&snes);
113:
115: /* setting the RHS function and the Jacobian's non-zero structutre */
116: SNESSetFunction(snes,global,FormFunction,&data);
117: SNESSetJacobian(snes,A,A,FormJacobian,&data);
119: /* set TSPVodeRHSFunction and TSPVodeRHSJacobian, so PETSc will pick up the
120: RHS function from SNES and compute the Jacobian by FD */
121: /*
122: TSSetRHSFunction(ts,TSPVodeSetRHSFunction,snes);
123: TSPVodeSetRHSJacobian(ts,0.0,global,&A,&A,&A_structure,snes);
124: TSSetRHSJacobian(ts,A,A,TSPVodeSetRHSJacobian,snes);
125: */
126:
127: TSSetRHSFunction(ts,RHSFunction,&data);
128: RHSJacobian(ts,0.0,global,&A,&A,&A_structure,&data);
129: TSSetRHSJacobian(ts,A,A,RHSJacobian,&data);
131: /* Use PVODE */
132: TSSetType(ts,TS_PVODE);
134: TSSetFromOptions(ts);
136: TSSetInitialTimeStep(ts,0.0,dt);
137: TSSetDuration(ts,time_steps,1);
138: TSSetSolution(ts,global);
141: /* Pick up a Petsc preconditioner */
142: /* one can always set method or preconditioner during the run time */
143: #if defined(PETSC_HAVE_PVODE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
144: TSPVodeGetPC(ts,&pc);
145: PCSetType(pc,PCJACOBI);
146: TSPVodeSetType(ts,PVODE_BDF);
147: /* TSPVodeSetMethodFromOptions(ts); */
148: #endif
150: TSSetUp(ts);
151: TSStep(ts,&steps,&ftime);
153: TSGetSolution(ts,&global);
154: PetscViewerASCIIOpen(PETSC_COMM_SELF,"out.m",&viewfile);
155: PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
156: VecView(global,viewfile);
158: #if defined(PETSC_HAVE_PVODE) && !defined(PETSC_USE_COMPLEX) && !defined(PETSC_USE_SINGLE)
159: /* extracts the PC from ts */
160: TSPVodeGetPC(ts,&pc);
161: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
162: TSView(ts,viewer);
163: PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
164: PCView(pc,viewer);
165: PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s Preconditioner,%s\n",
166: size,tsinfo,pcinfo);
167: PCDestroy(pc);
168: #endif
170: /* free the memories */
171: TSDestroy(ts);
172: VecDestroy(global);
173: if (A) {ierr= MatDestroy(A);}
174: PetscFinalize();
175: return 0;
176: }
178: /* -------------------------------------------------------------------*/
181: int Initial(Vec global,void *ctx)
182: {
183: Data *data = (Data*)ctx;
184: int m;
185: int row,col;
186: PetscReal x,y,dx,dy;
187: PetscScalar *localptr;
188: int i,mybase,myend,ierr,locsize;
190: /* make the local copies of parameters */
191: m = data->m;
192: dx = data->dx;
193: dy = data->dy;
195: /* determine starting point of each processor */
196: VecGetOwnershipRange(global,&mybase,&myend);
197: VecGetLocalSize(global,&locsize);
199: /* Initialize the array */
200: VecGetArray(global,&localptr);
202: for (i=0; i<locsize; i++) {
203: row = 1+(mybase+i)-((mybase+i)/m)*m;
204: col = (mybase+i)/m+1;
205: x = dx*row;
206: y = dy*col;
207: localptr[i] = f_ini(x,y);
208: }
209:
210: VecRestoreArray(global,&localptr);
211: return 0;
212: }
216: int Monitor(TS ts,int step,PetscReal time,Vec global,void *ctx)
217: {
218: VecScatter scatter;
219: IS from,to;
220: int i,n,*idx;
221: Vec tmp_vec;
222: int ierr;
223: PetscScalar *tmp;
225: /* Get the size of the vector */
226: VecGetSize(global,&n);
228: /* Set the index sets */
229: PetscMalloc(n*sizeof(int),&idx);
230: for(i=0; i<n; i++) idx[i]=i;
231:
232: /* Create local sequential vectors */
233: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
235: /* Create scatter context */
236: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
237: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
238: VecScatterCreate(global,from,tmp_vec,to,&scatter);
239: VecScatterBegin(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
240: VecScatterEnd(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
242: VecGetArray(tmp_vec,&tmp);
243: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u= %14.6e at the center \n",time,PetscRealPart(tmp[n/2]));
244: VecRestoreArray(tmp_vec,&tmp);
246: PetscFree(idx);
247: return 0;
248: }
253: int FormFunction(SNES snes,Vec globalin,Vec globalout,void *ptr)
254: {
255: Data *data = (Data*)ptr;
256: int m,n,mn;
257: PetscReal dx,dy;
258: PetscReal xc,xl,xr,yl,yr;
259: PetscReal a,epsilon;
260: PetscScalar *inptr,*outptr;
261: int i,j,len,ierr;
263: IS from,to;
264: int *idx;
265: VecScatter scatter;
266: Vec tmp_in,tmp_out;
268: m = data->m;
269: n = data->n;
270: mn = m*n;
271: dx = data->dx;
272: dy = data->dy;
273: a = data->a;
274: epsilon = data->epsilon;
276: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
277: xl = 0.5*a/dx+epsilon/(dx*dx);
278: xr = -0.5*a/dx+epsilon/(dx*dx);
279: yl = 0.5*a/dy+epsilon/(dy*dy);
280: yr = -0.5*a/dy+epsilon/(dy*dy);
281:
282: /* Get the length of parallel vector */
283: VecGetSize(globalin,&len);
285: /* Set the index sets */
286: PetscMalloc(len*sizeof(int),&idx);
287: for(i=0; i<len; i++) idx[i]=i;
288:
289: /* Create local sequential vectors */
290: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
291: VecDuplicate(tmp_in,&tmp_out);
293: /* Create scatter context */
294: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&from);
295: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&to);
296: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
297: VecScatterBegin(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
298:
300:
301: VecScatterEnd(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
302:
304: /*Extract income array */
305: VecGetArray(tmp_in,&inptr);
307: /* Extract outcome array*/
308: VecGetArray(tmp_out,&outptr);
310: outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
311: outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
312: for (i=1; i<m-1; i++) {
313: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
314: +yr*inptr[i+m];
315: }
317: for (j=1; j<n-1; j++) {
318: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
319: yl*inptr[j*m-m]+yr*inptr[j*m+m];
320: outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
321: yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
322: for(i=1; i<m-1; i++) {
323: outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
324: +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
325: }
326: }
328: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
329: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
330: for (i=1; i<m-1; i++) {
331: outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
332: +2*yl*inptr[mn-m+i-m];
333: }
335: VecRestoreArray(tmp_in,&inptr);
336: VecRestoreArray(tmp_out,&outptr);
338: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
339: VecScatterBegin(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter
340: );
341:
342: VecScatterEnd(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
343:
344:
345: /* Destroy idx aand scatter */
346: ISDestroy(from);
347: ISDestroy(to);
348: VecScatterDestroy(scatter);
350: PetscFree(idx);
351: return 0;
352: }
356: int FormJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
357: {
358: Data *data = (Data*)ptr;
359: Mat A = *AA;
360: PetscScalar v[1],one = 1.0;
361: int idx[1],i,j,row,ierr;
362: int m,n,mn;
364: m = data->m;
365: n = data->n;
366: mn = (data->m)*(data->n);
367:
368: for(i=0; i<mn; i++) {
369: idx[0] = i;
370: v[0] = one;
371: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
372: }
374: for(i=0; i<mn-m; i++) {
375: idx[0] = i+m;
376: v[0]= one;
377: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
378: }
380: for(i=m; i<mn; i++) {
381: idx[0] = i-m;
382: v[0] = one;
383: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
384: }
386: for(j=0; j<n; j++) {
387: for(i=0; i<m-1; i++) {
388: row = j*m+i;
389: idx[0] = j*m+i+1;
390: v[0] = one;
391: MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
392: }
393: }
395: for(j=0; j<n; j++) {
396: for(i=1; i<m; i++) {
397: row = j*m+i;
398: idx[0] = j*m+i-1;
399: v[0] = one;
400: MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
401: }
402: }
403:
404: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
405: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
408: *flag = SAME_NONZERO_PATTERN;
409: return 0;
410: }
414: int RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
415: {
416: Data *data = (Data*)ptr;
417: Mat A = *AA;
418: PetscScalar v[5];
419: int idx[5],i,j,row,ierr;
420: int m,n,mn;
421: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
423: m = data->m;
424: n = data->n;
425: mn = m*n;
426: dx = data->dx;
427: dy = data->dy;
428: a = data->a;
429: epsilon = data->epsilon;
431: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
432: xl = 0.5*a/dx+epsilon/(dx*dx);
433: xr = -0.5*a/dx+epsilon/(dx*dx);
434: yl = 0.5*a/dy+epsilon/(dy*dy);
435: yr = -0.5*a/dy+epsilon/(dy*dy);
437: row=0;
438: v[0] = xc; v[1]=xr; v[2]=yr;
439: idx[0]=0; idx[1]=2; idx[2]=m;
440: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
442: row=m-1;
443: v[0]=2.0*xl; v[1]=xc; v[2]=yr;
444: idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
445: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
447: for (i=1; i<m-1; i++) {
448: row=i;
449: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
450: idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
451: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
452: }
454: for (j=1; j<n-1; j++) {
455: row=j*m;
456: v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
457: idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
458: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
459:
460: row=j*m+m-1;
461: v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
462: 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;
463: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
465: for(i=1; i<m-1; i++) {
466: row=j*m+i;
467: v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
468: idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
469: idx[4]=j*m+i+m;
470: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
471: }
472: }
474: row=mn-m;
475: v[0] = xc; v[1]=xr; v[2]=2.0*yl;
476: idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
477: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
478:
479: row=mn-1;
480: v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
481: idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
482: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
484: for (i=1; i<m-1; i++) {
485: row=mn-m+i;
486: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
487: idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
488: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
489: }
492: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
493: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
496: *flag = SAME_NONZERO_PATTERN;
497: return 0;
498: }
502: int RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
503: {
505: SNES snes = PETSC_NULL;
507: FormFunction(snes,globalin,globalout,ctx);
510: return 0;
511: }