Actual source code: ex4.c
petsc-3.3-p6 2013-02-11
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;
31: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void *);
32: extern PetscErrorCode Initial(Vec,void *);
33: extern PetscErrorCode RHSFunction(TS,PetscReal,Vec,Vec,void*);
34: extern PetscErrorCode RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure *,void*);
35: extern PetscErrorCode PostStep(TS);
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);
88: #if defined(PETSC_HAVE_SUNDIALS)
89: TSSetType(ts,TSSUNDIALS);
90: #else
91: TSSetType(ts,TSEULER);
92: #endif
93: dt = 0.1;
94: ftime_original = data.tfinal = 1.0;
95: TSSetInitialTimeStep(ts,0.0,dt);
96: TSSetDuration(ts,time_steps,ftime_original);
97: TSSetSolution(ts,global);
99: /* set user provided RHSFunction and RHSJacobian */
100: TSSetRHSFunction(ts,PETSC_NULL,RHSFunction,&data);
101: MatCreate(PETSC_COMM_WORLD,&J);
102: MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
103: MatSetFromOptions(J);
104: MatSeqAIJSetPreallocation(J,5,PETSC_NULL);
105: MatMPIAIJSetPreallocation(J,5,PETSC_NULL,5,PETSC_NULL);
107: PetscOptionsHasName(PETSC_NULL,"-ts_fd",&flg);
108: if (!flg){
109: TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
110: } else {
111: TSGetSNES(ts,&snes);
112: PetscOptionsHasName(PETSC_NULL,"-fd_color",&fd_jacobian_coloring);
113: if (fd_jacobian_coloring){ /* Use finite differences with coloring */
114: /* Get data structure of J */
115: PetscBool pc_diagonal;
116: PetscOptionsHasName(PETSC_NULL,"-pc_diagonal",&pc_diagonal);
117: if (pc_diagonal){ /* the preconditioner of J is a diagonal matrix */
118: PetscInt rstart,rend,i;
119: PetscScalar zero=0.0;
120: MatGetOwnershipRange(J,&rstart,&rend);
121: for (i=rstart; i<rend; i++){
122: MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
123: }
124: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
125: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
126: } else {
127: /* Fill the structure using the expensive SNESDefaultComputeJacobian. Temporarily set up the TS so we can call this function */
128: TSSetType(ts,TSBEULER);
129: TSSetUp(ts);
130: SNESDefaultComputeJacobian(snes,x,&J,&J,&J_structure,ts);
131: }
133: /* create coloring context */
134: MatGetColoring(J,MATCOLORINGSL,&iscoloring);
135: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
136: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
137: MatFDColoringSetFromOptions(matfdcoloring);
138: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
139: ISColoringDestroy(&iscoloring);
140: } else { /* Use finite differences (slow) */
141: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobian,PETSC_NULL);
142: }
143: }
145: /* Pick up a Petsc preconditioner */
146: /* one can always set method or preconditioner during the run time */
147: TSGetSNES(ts,&snes);
148: SNESGetKSP(snes,&ksp);
149: KSPGetPC(ksp,&pc);
150: PCSetType(pc,PCJACOBI);
152: TSSetFromOptions(ts);
153: TSSetUp(ts);
154:
155: /* Test TSSetPostStep() */
156: PetscOptionsHasName(PETSC_NULL,"-test_PostStep",&flg);
157: if (flg){
158: TSSetPostStep(ts,PostStep);
159: }
161: PetscOptionsGetInt(PETSC_NULL,"-NOUT",&NOUT,PETSC_NULL);
162: for (iout=1; iout<=NOUT; iout++){
163: TSSetDuration(ts,time_steps,iout*ftime_original/NOUT);
164: TSSolve(ts,global,&ftime);
165: TSSetInitialTimeStep(ts,ftime,dt);
166: }
167: /* Interpolate solution at tfinal */
168: TSGetSolution(ts,&global);
169: TSInterpolate(ts,ftime_original,global);
171: PetscOptionsHasName(PETSC_NULL,"-matlab_view",&flg);
172: if (flg){ /* print solution into a MATLAB file */
173: PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
174: PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
175: VecView(global,viewfile);
176: PetscViewerDestroy(&viewfile);
177: }
179: /* display solver info for Sundials */
180: TSGetType(ts,&tstype);
181: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);
182: if (sundials){
183: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
184: TSView(ts,viewer);
185: PetscViewerDestroy(&viewer);
186: PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
187: PCView(pc,viewer);
188: PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",
189: size,tsinfo,pcinfo);
190: PetscViewerDestroy(&viewer);
191: }
193: /* free the memories */
194: TSDestroy(&ts);
195: VecDestroy(&global);
196: VecDestroy(&x);
197: ierr= MatDestroy(&J);
198: if (fd_jacobian_coloring){MatFDColoringDestroy(&matfdcoloring);}
199: PetscFinalize();
200: return 0;
201: }
203: /* -------------------------------------------------------------------*/
204: /* the initial function */
205: PetscReal f_ini(PetscReal x,PetscReal y)
206: {
207: PetscReal f;
209: f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
210: return f;
211: }
215: PetscErrorCode Initial(Vec global,void *ctx)
216: {
217: Data *data = (Data*)ctx;
218: PetscInt m,row,col;
219: PetscReal x,y,dx,dy;
220: PetscScalar *localptr;
221: PetscInt i,mybase,myend,locsize;
225: /* make the local copies of parameters */
226: m = data->m;
227: dx = data->dx;
228: dy = data->dy;
230: /* determine starting point of each processor */
231: VecGetOwnershipRange(global,&mybase,&myend);
232: VecGetLocalSize(global,&locsize);
234: /* Initialize the array */
235: VecGetArray(global,&localptr);
237: for (i=0; i<locsize; i++) {
238: row = 1+(mybase+i)-((mybase+i)/m)*m;
239: col = (mybase+i)/m+1;
240: x = dx*row;
241: y = dy*col;
242: localptr[i] = f_ini(x,y);
243: }
245: VecRestoreArray(global,&localptr);
246: return(0);
247: }
251: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
252: {
253: VecScatter scatter;
254: IS from,to;
255: PetscInt i,n,*idx,nsteps,maxsteps;
256: Vec tmp_vec;
258: PetscScalar *tmp;
259: PetscReal maxtime;
260: Data *data = (Data*)ctx;
261: PetscReal tfinal = data->tfinal;
264: if (time > tfinal) return(0);
266: TSGetTimeStepNumber(ts,&nsteps);
267: /* display output at selected time steps */
268: TSGetDuration(ts, &maxsteps, &maxtime);
269: if (nsteps % 10 != 0 && time < maxtime) return(0);
271: /* Get the size of the vector */
272: VecGetSize(global,&n);
274: /* Set the index sets */
275: PetscMalloc(n*sizeof(PetscInt),&idx);
276: for(i=0; i<n; i++) idx[i]=i;
278: /* Create local sequential vectors */
279: VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);
281: /* Create scatter context */
282: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
283: ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
284: VecScatterCreate(global,from,tmp_vec,to,&scatter);
285: VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
286: VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
288: VecGetArray(tmp_vec,&tmp);
289: PetscPrintf(PETSC_COMM_WORLD,"At t[%d] =%14.2e u= %14.2e at the center \n",nsteps,time,PetscRealPart(tmp[n/2]));
290: VecRestoreArray(tmp_vec,&tmp);
292: PetscFree(idx);
293: ISDestroy(&from);
294: ISDestroy(&to);
295: VecScatterDestroy(&scatter);
296: VecDestroy(&tmp_vec);
297: return(0);
298: }
302: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
303: {
304: Data *data = (Data*)ptr;
305: Mat A = *AA;
306: PetscScalar v[5];
307: PetscInt idx[5],i,j,row;
309: PetscInt m,n,mn;
310: PetscReal dx,dy,a,epsilon,xc,xl,xr,yl,yr;
313: m = data->m;
314: n = data->n;
315: mn = m*n;
316: dx = data->dx;
317: dy = data->dy;
318: a = data->a;
319: epsilon = data->epsilon;
321: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
322: xl = 0.5*a/dx+epsilon/(dx*dx);
323: xr = -0.5*a/dx+epsilon/(dx*dx);
324: yl = 0.5*a/dy+epsilon/(dy*dy);
325: yr = -0.5*a/dy+epsilon/(dy*dy);
327: row=0;
328: v[0] = xc; v[1]=xr; v[2]=yr;
329: idx[0]=0; idx[1]=2; idx[2]=m;
330: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
332: row=m-1;
333: v[0]=2.0*xl; v[1]=xc; v[2]=yr;
334: idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
335: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
337: for (i=1; i<m-1; i++) {
338: row=i;
339: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
340: idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
341: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
342: }
344: for (j=1; j<n-1; j++) {
345: row=j*m;
346: v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
347: idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
348: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
350: row=j*m+m-1;
351: v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
352: 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;
353: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
355: for(i=1; i<m-1; i++) {
356: row=j*m+i;
357: v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
358: idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
359: idx[4]=j*m+i+m;
360: MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
361: }
362: }
364: row=mn-m;
365: v[0] = xc; v[1]=xr; v[2]=2.0*yl;
366: idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
367: MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
369: row=mn-1;
370: v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
371: idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
372: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
374: for (i=1; i<m-1; i++) {
375: row=mn-m+i;
376: v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
377: idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
378: MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
379: }
381: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
382: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
384: /* *flag = SAME_NONZERO_PATTERN; */
385: *flag = DIFFERENT_NONZERO_PATTERN;
386: return(0);
387: }
389: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
392: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
393: {
394: Data *data = (Data*)ctx;
395: PetscInt m,n,mn;
396: PetscReal dx,dy;
397: PetscReal xc,xl,xr,yl,yr;
398: PetscReal a,epsilon;
399: PetscScalar *inptr,*outptr;
400: PetscInt i,j,len;
402: IS from,to;
403: PetscInt *idx;
404: VecScatter scatter;
405: Vec tmp_in,tmp_out;
408: m = data->m;
409: n = data->n;
410: mn = m*n;
411: dx = data->dx;
412: dy = data->dy;
413: a = data->a;
414: epsilon = data->epsilon;
416: xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
417: xl = 0.5*a/dx+epsilon/(dx*dx);
418: xr = -0.5*a/dx+epsilon/(dx*dx);
419: yl = 0.5*a/dy+epsilon/(dy*dy);
420: yr = -0.5*a/dy+epsilon/(dy*dy);
422: /* Get the length of parallel vector */
423: VecGetSize(globalin,&len);
425: /* Set the index sets */
426: PetscMalloc(len*sizeof(PetscInt),&idx);
427: for(i=0; i<len; i++) idx[i]=i;
429: /* Create local sequential vectors */
430: VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
431: VecDuplicate(tmp_in,&tmp_out);
433: /* Create scatter context */
434: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
435: ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
436: VecScatterCreate(globalin,from,tmp_in,to,&scatter);
437: VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
438: VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
439: VecScatterDestroy(&scatter);
441: /*Extract income array - include ghost points */
442: VecGetArray(tmp_in,&inptr);
444: /* Extract outcome array*/
445: VecGetArray(tmp_out,&outptr);
447: outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
448: outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
449: for (i=1; i<m-1; i++) {
450: outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
451: +yr*inptr[i+m];
452: }
454: for (j=1; j<n-1; j++) {
455: outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
456: yl*inptr[j*m-m]+yr*inptr[j*m+m];
457: outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
458: yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
459: for(i=1; i<m-1; i++) {
460: outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
461: +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
462: }
463: }
465: outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
466: outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
467: for (i=1; i<m-1; i++) {
468: outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
469: +2*yl*inptr[mn-m+i-m];
470: }
472: VecRestoreArray(tmp_in,&inptr);
473: VecRestoreArray(tmp_out,&outptr);
475: VecScatterCreate(tmp_out,from,globalout,to,&scatter);
476: VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
477: VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
479: /* Destroy idx aand scatter */
480: VecDestroy(&tmp_in);
481: VecDestroy(&tmp_out);
482: ISDestroy(&from);
483: ISDestroy(&to);
484: VecScatterDestroy(&scatter);
486: PetscFree(idx);
487: return(0);
488: }
492: PetscErrorCode PostStep(TS ts)
493: {
494: PetscErrorCode ierr;
495: PetscReal t;
498: TSGetTime(ts,&t);
499: PetscPrintf(PETSC_COMM_SELF," PostStep, t: %g\n",t);
500: return(0);
501: }