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:
10: This program tests the routine of computing the Jacobian by the
11: finite difference method as well as PETSc with SUNDIALS.
13: */
15: static char help[] = "Solve the convection-diffusion equation. \n\n";
17: #include petscsys.h
18: #include petscts.h
23: typedef struct
24: {
25: PetscInt m; /* the number of mesh points in x-direction */
26: PetscInt n; /* the number of mesh points in y-direction */
27: PetscReal dx; /* the grid space in x-direction */
28: PetscReal dy; /* the grid space in y-direction */
29: PetscReal a; /* the convection coefficient */
30: PetscReal epsilon; /* the diffusion coefficient */
31: } Data;
33: /* two temporal functions */
39: /* the initial function */
40: PetscReal f_ini(PetscReal x,PetscReal y)
41: {
42: PetscReal f;
43: f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
44: return f;
45: }
47: /* #undef PETSC_HAVE_SUNDIALS */
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: {
60: PetscInt time_steps = 100,steps;
61: PetscMPIInt size;
62: Vec global;
63: PetscReal dt,ftime;
64: TS ts;
65: PetscViewer viewfile;
66: MatStructure A_structure;
67: Mat A = 0;
68: TSProblemType tsproblem = TS_NONLINEAR; /* Need to be TS_NONLINEAR */
69: SNES snes;
70: Vec x;
71: Data data;
72: PetscInt mn;
73: #if defined(PETSC_HAVE_SUNDIALS)
74: PC pc;
75: PetscViewer viewer;
76: char pcinfo[120],tsinfo[120];
77: #endif
79: PetscInitialize(&argc,&argv,(char*)0,help);
80: MPI_Comm_size(PETSC_COMM_WORLD,&size);
81:
82: /* set Data */
83: data.m = 9;
84: data.n = 9;
85: data.a = 1.0;
86: data.epsilon = 0.1;
87: data.dx = 1.0/(data.m+1.0);
88: data.dy = 1.0/(data.n+1.0);
89: mn = (data.m)*(data.n);
91: PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
92:
93: /* set initial conditions */
94: VecCreate(PETSC_COMM_WORLD,&global);
95: VecSetSizes(global,PETSC_DECIDE,mn);
96: VecSetFromOptions(global);
97: Initial(global,&data);
98: VecDuplicate(global,&x);
99:
100: /* make timestep context */
101: TSCreate(PETSC_COMM_WORLD,&ts);
102: TSSetProblemType(ts,tsproblem);
103: TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);
105: dt = 0.1;
107: /*
108: The user provides the RHS and Jacobian
109: */
110: MatCreate(PETSC_COMM_WORLD,&A);
111: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
112: MatSetFromOptions(A);
114: /* Create SNES context */
115: SNESCreate(PETSC_COMM_WORLD,&snes);
117: /* setting the RHS function and the Jacobian's non-zero structutre */
118: SNESSetFunction(snes,global,FormFunction,&data);
119: SNESSetJacobian(snes,A,A,FormJacobian,&data);
121: /* set TSSundialsRHSFunction and TSSundialsRHSJacobian, so PETSc will pick up the
122: RHS function from SNES and compute the Jacobian by FD */
123: /*
124: TSSetRHSFunction(ts,TSSundialsSetRHSFunction,snes);
125: TSSundialsSetRHSJacobian(ts,0.0,global,&A,&A,&A_structure,snes);
126: TSSetRHSJacobian(ts,A,A,TSSundialsSetRHSJacobian,snes);
127: */
128:
129: TSSetRHSFunction(ts,RHSFunction,&data);
130: RHSJacobian(ts,0.0,global,&A,&A,&A_structure,&data);
131: TSSetRHSJacobian(ts,A,A,RHSJacobian,&data);
133: /* Use SUNDIALS */
134: #if defined(PETSC_HAVE_SUNDIALS)
135: TSSetType(ts,TS_SUNDIALS);
136: #else
137: TSSetType(ts,TS_EULER);
138: #endif
140: TSSetFromOptions(ts);
142: TSSetInitialTimeStep(ts,0.0,dt);
143: TSSetDuration(ts,time_steps,1);
144: TSSetSolution(ts,global);
147: /* Pick up a Petsc preconditioner */
148: /* one can always set method or preconditioner during the run time */
149: #if defined(PETSC_HAVE_SUNDIALS)
150: TSSundialsGetPC(ts,&pc);
151: PCSetType(pc,PCJACOBI);
152: TSSundialsSetType(ts,SUNDIALS_BDF);
153: /* TSSundialsSetMethodFromOptions(ts); */
154: #endif
156: TSSetUp(ts);
157: TSStep(ts,&steps,&ftime);
159: TSGetSolution(ts,&global);
160: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
161: PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
162: VecView(global,viewfile);
163: PetscViewerDestroy(viewfile);
165: #if defined(PETSC_HAVE_SUNDIALS)
166: /* extracts the PC from ts */
167: TSSundialsGetPC(ts,&pc);
168: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
169: TSView(ts,viewer);
170: PetscViewerDestroy(viewer);
171: PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
172: PCView(pc,viewer);
173: PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s Preconditioner,%s\n",
174: size,tsinfo,pcinfo);
175: PetscViewerDestroy(viewer);
176: #endif
178: /* free the memories */
179: TSDestroy(ts);
180: VecDestroy(global);
181: VecDestroy(x);
182: if (A) {ierr= MatDestroy(A);}
183: SNESDestroy(snes);
184: PetscFinalize();
185: return 0;
186: }
188: /* -------------------------------------------------------------------*/
191: PetscErrorCode Initial(Vec global,void *ctx)
192: {
193: Data *data = (Data*)ctx;
194: PetscInt m,row,col;
195: PetscReal x,y,dx,dy;
196: PetscScalar *localptr;
197: PetscInt i,mybase,myend,locsize;
200: /* make the local copies of parameters */
201: m = data->m;
202: dx = data->dx;
203: dy = data->dy;
205: /* determine starting point of each processor */
206: VecGetOwnershipRange(global,&mybase,&myend);
207: VecGetLocalSize(global,&locsize);
209: /* Initialize the array */
210: VecGetArray(global,&localptr);
212: for (i=0; i<locsize; i++) {
213: row = 1+(mybase+i)-((mybase+i)/m)*m;
214: col = (mybase+i)/m+1;
215: x = dx*row;
216: y = dy*col;
217: localptr[i] = f_ini(x,y);
218: }
219:
220: VecRestoreArray(global,&localptr);
221: return 0;
222: }
226: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
227: {
228: VecScatter scatter;
229: IS from,to;
230: PetscInt i,n,*idx;
231: Vec tmp_vec;
233: PetscScalar *tmp;
235: /* Get the size of the vector */
236: VecGetSize(global,&n);
238: /* Set the index sets */
239: PetscMalloc(n*sizeof(PetscInt),&idx);
240: for(i=0; i<n; i++) idx[i]=i;
241:
242: /* Create local sequential vectors */
243: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
245: /* Create scatter context */
246: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
247: ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
248: VecScatterCreate(global,from,tmp_vec,to,&scatter);
249: VecScatterBegin(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
250: VecScatterEnd(global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD,scatter);
252: VecGetArray(tmp_vec,&tmp);
253: PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u= %14.6e at the center \n",time,PetscRealPart(tmp[n/2]));
254: VecRestoreArray(tmp_vec,&tmp);
256: PetscFree(idx);
257: ISDestroy(from);
258: ISDestroy(to);
259: VecScatterDestroy(scatter);
260: VecDestroy(tmp_vec);
261: return 0;
262: }
267: PetscErrorCode FormFunction(SNES snes,Vec globalin,Vec globalout,void *ptr)
268: {
269: Data *data = (Data*)ptr;
270: PetscInt m,n,mn;
271: PetscReal dx,dy;
272: PetscReal xc,xl,xr,yl,yr;
273: PetscReal a,epsilon;
274: PetscScalar *inptr,*outptr;
275: PetscInt i,j,len;
277: IS from,to;
278: PetscInt *idx;
279: VecScatter scatter;
280: Vec tmp_in,tmp_out;
282: m = data->m;
283: n = data->n;
284: mn = m*n;
285: dx = data->dx;
286: dy = data->dy;
287: a = data->a;
288: epsilon = data->epsilon;
290: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
291: xl = 0.5*a/dx+epsilon/(dx*dx);
292: xr = -0.5*a/dx+epsilon/(dx*dx);
293: yl = 0.5*a/dy+epsilon/(dy*dy);
294: yr = -0.5*a/dy+epsilon/(dy*dy);
295:
296: /* Get the length of parallel vector */
297: VecGetSize(globalin,&len);
299: /* Set the index sets */
300: PetscMalloc(len*sizeof(PetscInt),&idx);
301: for(i=0; i<len; i++) idx[i]=i;
302:
303: /* Create local sequential vectors */
304: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
305: VecDuplicate(tmp_in,&tmp_out);
307: /* Create scatter context */
308: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&from);
309: ISCreateGeneral(PETSC_COMM_SELF,len,idx,&to);
310: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
311: VecScatterBegin(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
312: VecScatterEnd(globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD,scatter);
313: VecScatterDestroy(scatter);
315: /*Extract income array */
316: VecGetArray(tmp_in,&inptr);
318: /* Extract outcome array*/
319: VecGetArray(tmp_out,&outptr);
321: outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
322: outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
323: for (i=1; i<m-1; i++) {
324: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
325: +yr*inptr[i+m];
326: }
328: for (j=1; j<n-1; j++) {
329: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
330: yl*inptr[j*m-m]+yr*inptr[j*m+m];
331: outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
332: yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
333: for(i=1; i<m-1; i++) {
334: outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
335: +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
336: }
337: }
339: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
340: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
341: for (i=1; i<m-1; i++) {
342: outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
343: +2*yl*inptr[mn-m+i-m];
344: }
346: VecRestoreArray(tmp_in,&inptr);
347: VecRestoreArray(tmp_out,&outptr);
349: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
350: VecScatterBegin(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
351: VecScatterEnd(tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD,scatter);
352:
353: /* Destroy idx aand scatter */
354: VecDestroy(tmp_in);
355: VecDestroy(tmp_out);
356: ISDestroy(from);
357: ISDestroy(to);
358: VecScatterDestroy(scatter);
360: PetscFree(idx);
361: return 0;
362: }
366: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
367: {
368: Data *data = (Data*)ptr;
369: Mat A = *AA;
370: PetscScalar v[1],one = 1.0;
371: PetscInt idx[1],i,j,row;
373: PetscInt m,n,mn;
375: m = data->m;
376: n = data->n;
377: mn = (data->m)*(data->n);
378:
379: for(i=0; i<mn; i++) {
380: idx[0] = i;
381: v[0] = one;
382: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
383: }
385: for(i=0; i<mn-m; i++) {
386: idx[0] = i+m;
387: v[0]= one;
388: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
389: }
391: for(i=m; i<mn; i++) {
392: idx[0] = i-m;
393: v[0] = one;
394: MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
395: }
397: for(j=0; j<n; j++) {
398: for(i=0; i<m-1; i++) {
399: row = j*m+i;
400: idx[0] = j*m+i+1;
401: v[0] = one;
402: MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
403: }
404: }
406: for(j=0; j<n; j++) {
407: for(i=1; i<m; i++) {
408: row = j*m+i;
409: idx[0] = j*m+i-1;
410: v[0] = one;
411: MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
412: }
413: }
414:
415: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
416: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
419: *flag = SAME_NONZERO_PATTERN;
420: return 0;
421: }
425: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
426: {
427: Data *data = (Data*)ptr;
428: Mat A = *AA;
429: PetscScalar v[5];
430: PetscInt idx[5],i,j,row;
432: PetscInt m,n,mn;
433: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
435: m = data->m;
436: n = data->n;
437: mn = m*n;
438: dx = data->dx;
439: dy = data->dy;
440: a = data->a;
441: epsilon = data->epsilon;
443: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
444: xl = 0.5*a/dx+epsilon/(dx*dx);
445: xr = -0.5*a/dx+epsilon/(dx*dx);
446: yl = 0.5*a/dy+epsilon/(dy*dy);
447: yr = -0.5*a/dy+epsilon/(dy*dy);
449: row=0;
450: v[0] = xc; v[1]=xr; v[2]=yr;
451: idx[0]=0; idx[1]=2; idx[2]=m;
452: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
454: row=m-1;
455: v[0]=2.0*xl; v[1]=xc; v[2]=yr;
456: idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
457: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
459: for (i=1; i<m-1; i++) {
460: row=i;
461: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
462: idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
463: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
464: }
466: for (j=1; j<n-1; j++) {
467: row=j*m;
468: v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
469: idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
470: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
471:
472: row=j*m+m-1;
473: v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
474: 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;
475: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
477: for(i=1; i<m-1; i++) {
478: row=j*m+i;
479: v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
480: idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
481: idx[4]=j*m+i+m;
482: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
483: }
484: }
486: row=mn-m;
487: v[0] = xc; v[1]=xr; v[2]=2.0*yl;
488: idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
489: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
490:
491: row=mn-1;
492: v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
493: idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
494: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
496: for (i=1; i<m-1; i++) {
497: row=mn-m+i;
498: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
499: idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
500: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
501: }
504: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
505: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
508: /* *flag = SAME_NONZERO_PATTERN; */
509: *flag = DIFFERENT_NONZERO_PATTERN;
510: return 0;
511: }
515: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
516: {
518: SNES snes = PETSC_NULL;
520: FormFunction(snes,globalin,globalout,ctx);
523: return 0;
524: }