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: }