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;


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

 64:   PetscInitialize(&argc,&argv,(char*)0,help);
 65:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
 66: 
 67:   /* set data */
 68:   data.m = 9;
 69:   data.n = 9;
 70:   data.a = 1.0;
 71:   data.epsilon = 0.1;
 72:   data.dx      = 1.0/(data.m+1.0);
 73:   data.dy      = 1.0/(data.n+1.0);
 74:   mn = (data.m)*(data.n);
 75:   PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
 76: 
 77:   /* set initial conditions */
 78:   VecCreate(PETSC_COMM_WORLD,&global);
 79:   VecSetSizes(global,PETSC_DECIDE,mn);
 80:   VecSetFromOptions(global);
 81:   Initial(global,&data);
 82:   VecDuplicate(global,&x);

 84:   /* create timestep context */
 85:   TSCreate(PETSC_COMM_WORLD,&ts);
 86:   TSSetProblemType(ts,TS_NONLINEAR); /* Need to be TS_NONLINEAR for Sundials */
 87:   TSMonitorSet(ts,Monitor,PETSC_NULL,PETSC_NULL);

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

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

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

140:   /* Use SUNDIALS */
141: #if defined(PETSC_HAVE_SUNDIALS)
142:   TSSetType(ts,TSSUNDIALS);
143: #else
144:   TSSetType(ts,TSEULER);
145: #endif
146:   dt   = 0.1;
147:   TSSetInitialTimeStep(ts,0.0,dt);
148:   TSSetDuration(ts,time_steps,1);
149:   TSSetSolution(ts,global);
150: 
151:   /* Test TSSetPostStep() */
152:   PetscOptionsHasName(PETSC_NULL,"-test_PostStep",&flg);
153:   if (flg){
154:     TSSetPostStep(ts,PostStep);
155:   }

157:   /* Pick up a Petsc preconditioner */
158:   /* one can always set method or preconditioner during the run time */
159: #if defined(PETSC_HAVE_SUNDIALS)
160:   TSSundialsGetPC(ts,&pc);
161:   PCSetType(pc,PCJACOBI);
162: #endif
163:   TSSetFromOptions(ts);
164:   TSSetUp(ts);

166:   PetscOptionsGetInt(PETSC_NULL,"-NOUT",&NOUT,PETSC_NULL);
167:   for (iout=1; iout<=NOUT; iout++){
168:     TSSetDuration(ts,time_steps,iout*1.0/NOUT);
169:     TSStep(ts,&steps,&ftime);
170:     TSSetInitialTimeStep(ts,ftime,dt);
171:   }

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,TSSUNDIALS,&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,maxsteps;
262:   Vec            tmp_vec;
264:   PetscScalar    *tmp;
265:   PetscReal      maxtime;
266: 
268:   TSGetTimeStepNumber(ts,&nsteps);
269:   /* display output at selected time steps */
270:   TSGetDuration(ts, &maxsteps, &maxtime);
271:   if (nsteps % 10 != 0 && time < maxtime) return(0);

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

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

283:   /* Create scatter context */
284:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
285:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
286:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
287:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
288:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

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

294:   PetscFree(idx);
295:   ISDestroy(from);
296:   ISDestroy(to);
297:   VecScatterDestroy(scatter);
298:   VecDestroy(tmp_vec);
299:   return(0);
300: }

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

315:   m = data->m;
316:   n = data->n;
317:   mn = m*n;
318:   dx = data->dx;
319:   dy = data->dy;
320:   a = data->a;
321:   epsilon = data->epsilon;

323:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
324:   xl = 0.5*a/dx+epsilon/(dx*dx);
325:   xr = -0.5*a/dx+epsilon/(dx*dx);
326:   yl = 0.5*a/dy+epsilon/(dy*dy);
327:   yr = -0.5*a/dy+epsilon/(dy*dy);

329:   row=0;
330:   v[0] = xc; v[1]=xr; v[2]=yr;
331:   idx[0]=0; idx[1]=2; idx[2]=m;
332:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

334:   row=m-1;
335:   v[0]=2.0*xl; v[1]=xc; v[2]=yr;
336:   idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
337:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

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

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

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

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

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

383:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
384:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

386:   /* *flag = SAME_NONZERO_PATTERN; */
387:   *flag = DIFFERENT_NONZERO_PATTERN;
388:   return(0);
389: }

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

410:   m       = data->m;
411:   n       = data->n;
412:   mn      = m*n;
413:   dx      = data->dx;
414:   dy      = data->dy;
415:   a       = data->a;
416:   epsilon = data->epsilon;

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

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

435:   /* Create scatter context */
436:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,&from);
437:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,&to);
438:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
439:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
440:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
441:   VecScatterDestroy(scatter);

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

446:   /* Extract outcome array*/
447:   VecGetArray(tmp_out,&outptr);

449:   outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
450:   outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
451:   for (i=1; i<m-1; i++) {
452:     outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
453:       +yr*inptr[i+m];
454:   }

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

467:   outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
468:   outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
469:   for (i=1; i<m-1; i++) {
470:     outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
471:       +2*yl*inptr[mn-m+i-m];
472:   }

474:   VecRestoreArray(tmp_in,&inptr);
475:   VecRestoreArray(tmp_out,&outptr);

477:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
478:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
479:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
480: 
481:   /* Destroy idx aand scatter */
482:   VecDestroy(tmp_in);
483:   VecDestroy(tmp_out);
484:   ISDestroy(from);
485:   ISDestroy(to);
486:   VecScatterDestroy(scatter);

488:   PetscFree(idx);
489:   return(0);
490: }

494: PetscErrorCode PostStep(TS ts)
495: {
496:   PetscErrorCode  ierr;
497:   PetscReal       t;
498: 
500:   TSGetTime(ts,&t);
501:   PetscPrintf(PETSC_COMM_SELF,"  PostStep, t: %g\n",t);
502:   return(0);
503: }