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


 39: int main(int argc,char **argv)
 40: {
 42:   PetscInt       time_steps = 100,steps;
 43:   PetscMPIInt    size;
 44:   Vec            global;
 45:   PetscReal      dt,ftime;
 46:   TS             ts;
 47:   PetscViewer         viewfile;
 48:   MatStructure   J_structure;
 49:   Mat            J = 0;
 50:   Vec                  x;
 51:   Data                 data;
 52:   PetscInt          mn;
 53:   PetscTruth     flg;
 54:   ISColoring     iscoloring;
 55:   MatFDColoring  matfdcoloring = 0;
 56:   PetscTruth     fd_jacobian_coloring = PETSC_FALSE;
 57: #if defined(PETSC_HAVE_SUNDIALS)
 58:   PC                 pc;
 59:   PetscViewer    viewer;
 60:   char           pcinfo[120],tsinfo[120];
 61:   const TSType   tstype;
 62:   PetscTruth     sundials;
 63: #endif

 65:   PetscInitialize(&argc,&argv,(char*)0,help);
 66:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 67: 
 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);
 77: 
 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:   TSSetProblemType(ts,TS_NONLINEAR); /* Need to be TS_NONLINEAR for Sundials */
 88:   TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);

 90:   /* set user provided RHSFunction and RHSJacobian */
 91:   TSSetRHSFunction(ts,RHSFunction,&data);
 92:   MatCreate(PETSC_COMM_WORLD,&J);
 93:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
 94:   MatSetFromOptions(J);
 95:   MatSeqAIJSetPreallocation(J,5,PETSC_NULL);
 96:   MatMPIAIJSetPreallocation(J,5,PETSC_NULL,5,PETSC_NULL);

 98:   PetscOptionsHasName(PETSC_NULL,"-ts_fd",&flg);
 99:   if (!flg){
100:     /* RHSJacobian(ts,0.0,global,&J,&J,&J_structure,&data); */
101:     TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
102:   } else {
103:     PetscOptionsHasName(PETSC_NULL,"-fd_color",&fd_jacobian_coloring);
104:     if (fd_jacobian_coloring){ /* Use finite differences with coloring */
105:       /* Get data structure of J */
106:       PetscTruth pc_diagonal;
107:       PetscOptionsHasName(PETSC_NULL,"-pc_diagonal",&pc_diagonal);
108:       if (pc_diagonal){ /* the preconditioner of J is a diagonal matrix */
109:         PetscInt rstart,rend,i;
110:         PetscScalar  zero=0.0;
111:         MatGetOwnershipRange(J,&rstart,&rend);
112:         for (i=rstart; i<rend; i++){
113:           MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
114:         }
115:         MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
116:         MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
117:       } else { /* the preconditioner of J has same data structure as J */
118:         TSDefaultComputeJacobian(ts,0.0,global,&J,&J,&J_structure,&data);
119:       }
120: 
121:       /* create coloring context */
122:       MatGetColoring(J,MATCOLORING_SL,&iscoloring);
123:       MatFDColoringCreate(J,iscoloring,&matfdcoloring);
124:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))RHSFunction,&data);
125:       MatFDColoringSetFromOptions(matfdcoloring);
126:       TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobianColor,matfdcoloring);
127:       ISColoringDestroy(iscoloring);

129:       PetscTruth view_J;
130:       PetscOptionsHasName(PETSC_NULL,"-view_J",&view_J);
131:       if (view_J){
132:         PetscPrintf(PETSC_COMM_SELF,"J computed from TSDefaultComputeJacobianColor():\n");
133:         TSDefaultComputeJacobianColor(ts,0.0,global,&J,&J,&J_structure,matfdcoloring);
134:         MatView(J,PETSC_VIEWER_STDOUT_WORLD);
135:       }
136:     } else { /* Use finite differences (slow) */
137:       TSSetRHSJacobian(ts,J,J,TSDefaultComputeJacobian,&data);
138:     }
139:   }

141:   /* Use SUNDIALS */
142: #if defined(PETSC_HAVE_SUNDIALS)
143:   TSSetType(ts,TS_SUNDIALS);
144: #else
145:   TSSetType(ts,TS_EULER);
146: #endif
147:   dt   = 0.1;
148:   TSSetInitialTimeStep(ts,0.0,dt);
149:   TSSetDuration(ts,time_steps,1);
150:   TSSetSolution(ts,global);
151: 
152:   /* Test TSSetPostStep() and TSSetPostUpdate() */
153:   PetscOptionsHasName(PETSC_NULL,"-test_PostStep",&flg);
154:   if (flg){
155:     TSSetPostStep(ts,PostStep);
156:   }
157:   PetscOptionsHasName(PETSC_NULL,"-test_PostUpdate",&flg);
158:   if (flg){
159:     TSSetPostUpdate(ts,PostUpdate);
160:   }

162:   /* Pick up a Petsc preconditioner */
163:   /* one can always set method or preconditioner during the run time */
164: #if defined(PETSC_HAVE_SUNDIALS)
165:   TSSundialsGetPC(ts,&pc);
166:   PCSetType(pc,PCJACOBI);
167: #endif
168:   TSSetFromOptions(ts);

170:   TSSetUp(ts);
171:   TSStep(ts,&steps,&ftime);

173:   PetscOptionsHasName(PETSC_NULL,"-matlab_view",&flg);
174:   if (flg){ /* print solution into a MATLAB file */
175:     TSGetSolution(ts,&global);
176:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
177:     PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
178:     VecView(global,viewfile);
179:     PetscViewerDestroy(viewfile);
180:   }

182: #if defined(PETSC_HAVE_SUNDIALS)
183:   /* extracts the PC  from ts */
184:   TSSundialsGetPC(ts,&pc);
185:   TSGetType(ts,&tstype);
186:   PetscTypeCompare((PetscObject)ts,TS_SUNDIALS,&sundials);
187:   if (sundials){
188:     PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
189:     TSView(ts,viewer);
190:     PetscViewerDestroy(viewer);
191:     PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
192:     PCView(pc,viewer);
193:     PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",
194:                      size,tsinfo,pcinfo);
195:     PetscViewerDestroy(viewer);
196:   }
197: #endif

199:   /* free the memories */
200:   TSDestroy(ts);
201:   VecDestroy(global);
202:   VecDestroy(x);
203:   ierr= MatDestroy(J);
204:   if (fd_jacobian_coloring){MatFDColoringDestroy(matfdcoloring);}
205:   PetscFinalize();
206:   return 0;
207: }

209: /* -------------------------------------------------------------------*/
210: /* the initial function */
211: PetscReal f_ini(PetscReal x,PetscReal y)
212: {
213:   PetscReal f;

215:   f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
216:   return f;
217: }

221: PetscErrorCode Initial(Vec global,void *ctx)
222: {
223:   Data           *data = (Data*)ctx;
224:   PetscInt       m,row,col;
225:   PetscReal      x,y,dx,dy;
226:   PetscScalar    *localptr;
227:   PetscInt       i,mybase,myend,locsize;

231:   /* make the local  copies of parameters */
232:   m = data->m;
233:   dx = data->dx;
234:   dy = data->dy;

236:   /* determine starting point of each processor */
237:   VecGetOwnershipRange(global,&mybase,&myend);
238:   VecGetLocalSize(global,&locsize);

240:   /* Initialize the array */
241:   VecGetArray(global,&localptr);

243:   for (i=0; i<locsize; i++) {
244:     row = 1+(mybase+i)-((mybase+i)/m)*m;
245:     col = (mybase+i)/m+1;
246:     x = dx*row;
247:     y = dy*col;
248:     localptr[i] = f_ini(x,y);
249:   }
250: 
251:   VecRestoreArray(global,&localptr);
252:   return(0);
253: }

257: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
258: {
259:   VecScatter     scatter;
260:   IS             from,to;
261:   PetscInt       i,n,*idx,nsteps;
262:   Vec            tmp_vec;
264:   PetscScalar    *tmp;
265: 
267:   TSGetTimeStepNumber(ts,&nsteps);

269:   /* Get the size of the vector */
270:   VecGetSize(global,&n);

272:   /* Set the index sets */
273:   PetscMalloc(n*sizeof(PetscInt),&idx);
274:   for(i=0; i<n; i++) idx[i]=i;
275: 
276:   /* Create local sequential vectors */
277:   VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);

279:   /* Create scatter context */
280:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
281:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
282:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
283:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
284:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

286:   VecGetArray(tmp_vec,&tmp);
287:   PetscPrintf(PETSC_COMM_WORLD,"At t[%d] =%14.2e u= %14.2e at the center \n",nsteps,time,PetscRealPart(tmp[n/2]));
288:   VecRestoreArray(tmp_vec,&tmp);

290:   PetscFree(idx);
291:   ISDestroy(from);
292:   ISDestroy(to);
293:   VecScatterDestroy(scatter);
294:   VecDestroy(tmp_vec);
295:   return(0);
296: }

300: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
301: {
302:   Data           *data = (Data*)ptr;
303:   Mat            A = *AA;
304:   PetscScalar    v[5];
305:   PetscInt       idx[5],i,j,row;
307:   PetscInt       m,n,mn;
308:   PetscReal      dx,dy,a,epsilon,xc,xl,xr,yl,yr;

311:   m = data->m;
312:   n = data->n;
313:   mn = m*n;
314:   dx = data->dx;
315:   dy = data->dy;
316:   a = data->a;
317:   epsilon = data->epsilon;

319:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
320:   xl = 0.5*a/dx+epsilon/(dx*dx);
321:   xr = -0.5*a/dx+epsilon/(dx*dx);
322:   yl = 0.5*a/dy+epsilon/(dy*dy);
323:   yr = -0.5*a/dy+epsilon/(dy*dy);

325:   row=0;
326:   v[0] = xc; v[1]=xr; v[2]=yr;
327:   idx[0]=0; idx[1]=2; idx[2]=m;
328:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

330:   row=m-1;
331:   v[0]=2.0*xl; v[1]=xc; v[2]=yr;
332:   idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
333:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

335:   for (i=1; i<m-1; i++) {
336:     row=i;
337:     v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
338:     idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
339:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
340:   }

342:   for (j=1; j<n-1; j++) {
343:     row=j*m;
344:     v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
345:     idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
346:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
347: 
348:     row=j*m+m-1;
349:     v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
350:     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;
351:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);

353:     for(i=1; i<m-1; i++) {
354:       row=j*m+i;
355:       v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
356:       idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
357:       idx[4]=j*m+i+m;
358:       MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
359:     }
360:   }

362:   row=mn-m;
363:   v[0] = xc; v[1]=xr; v[2]=2.0*yl;
364:   idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
365:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
366: 
367:   row=mn-1;
368:   v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
369:   idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
370:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);

372:   for (i=1; i<m-1; i++) {
373:     row=mn-m+i;
374:     v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
375:     idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
376:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
377:   }

379:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
380:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

382:   /* *flag = SAME_NONZERO_PATTERN; */
383:   *flag = DIFFERENT_NONZERO_PATTERN;
384:   return(0);
385: }

387: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
390: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
391: {
392:   Data           *data = (Data*)ctx;
393:   PetscInt       m,n,mn;
394:   PetscReal      dx,dy;
395:   PetscReal      xc,xl,xr,yl,yr;
396:   PetscReal      a,epsilon;
397:   PetscScalar    *inptr,*outptr;
398:   PetscInt       i,j,len;
400:   IS             from,to;
401:   PetscInt       *idx;
402:   VecScatter     scatter;
403:   Vec            tmp_in,tmp_out;

406:   m       = data->m;
407:   n       = data->n;
408:   mn      = m*n;
409:   dx      = data->dx;
410:   dy      = data->dy;
411:   a       = data->a;
412:   epsilon = data->epsilon;

414:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
415:   xl = 0.5*a/dx+epsilon/(dx*dx);
416:   xr = -0.5*a/dx+epsilon/(dx*dx);
417:   yl = 0.5*a/dy+epsilon/(dy*dy);
418:   yr = -0.5*a/dy+epsilon/(dy*dy);
419: 
420:   /* Get the length of parallel vector */
421:   VecGetSize(globalin,&len);

423:   /* Set the index sets */
424:   PetscMalloc(len*sizeof(PetscInt),&idx);
425:   for(i=0; i<len; i++) idx[i]=i;
426: 
427:   /* Create local sequential vectors */
428:   VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
429:   VecDuplicate(tmp_in,&tmp_out);

431:   /* Create scatter context */
432:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,&from);
433:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,&to);
434:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
435:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
436:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
437:   VecScatterDestroy(scatter);

439:   /*Extract income array - include ghost points */
440:   VecGetArray(tmp_in,&inptr);

442:   /* Extract outcome array*/
443:   VecGetArray(tmp_out,&outptr);

445:   outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
446:   outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
447:   for (i=1; i<m-1; i++) {
448:     outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
449:       +yr*inptr[i+m];
450:   }

452:   for (j=1; j<n-1; j++) {
453:     outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
454:       yl*inptr[j*m-m]+yr*inptr[j*m+m];
455:     outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
456:       yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
457:     for(i=1; i<m-1; i++) {
458:       outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
459:         +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
460:     }
461:   }

463:   outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
464:   outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
465:   for (i=1; i<m-1; i++) {
466:     outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
467:       +2*yl*inptr[mn-m+i-m];
468:   }

470:   VecRestoreArray(tmp_in,&inptr);
471:   VecRestoreArray(tmp_out,&outptr);

473:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
474:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
475:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
476: 
477:   /* Destroy idx aand scatter */
478:   VecDestroy(tmp_in);
479:   VecDestroy(tmp_out);
480:   ISDestroy(from);
481:   ISDestroy(to);
482:   VecScatterDestroy(scatter);

484:   PetscFree(idx);
485:   return(0);
486: }

490: PetscErrorCode PostStep(TS ts)
491: {
492:   PetscErrorCode  ierr;
493:   PetscReal       t;
494: 
496:   TSGetTime(ts,&t);
497:   PetscPrintf(PETSC_COMM_SELF,"  PostStep, t: %g\n",t);
498:   return(0);
499: }

503: PetscErrorCode PostUpdate(TS ts,PetscReal t, PetscReal *dt)
504: {
505:   PetscErrorCode  ierr;
507:   PetscPrintf(PETSC_COMM_SELF,"  PostUpdate, t: %g\n",t);
508:   return(0);
509: }