Actual source code: ex4.c

  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;


 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);

 89:   /* set user provided RHSFunction and RHSJacobian */
 90:   TSSetRHSFunction(ts,PETSC_NULL,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:       PetscBool  pc_diagonal;
106:       PetscBool  view_J;
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:       }

121:       /* create coloring context */
122:       MatGetColoring(J,MATCOLORINGSL,&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:       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:   ftime_original = data.tfinal = 1.0;
148:   TSSetInitialTimeStep(ts,0.0,dt);
149:   TSSetDuration(ts,time_steps,ftime_original);
150:   TSSetSolution(ts,global);


153:   /* Pick up a Petsc preconditioner */
154:   /* one can always set method or preconditioner during the run time */
155:   TSGetSNES(ts,&snes);
156:   SNESGetKSP(snes,&ksp);
157:   KSPGetPC(ksp,&pc);
158:   PCSetType(pc,PCJACOBI);

160:   TSSetFromOptions(ts);
161:   TSSetUp(ts);
162: 
163:   /* Test TSSetPostStep() */
164:   PetscOptionsHasName(PETSC_NULL,"-test_PostStep",&flg);
165:   if (flg){
166:     TSSetPostStep(ts,PostStep);
167:   }

169:   PetscOptionsGetInt(PETSC_NULL,"-NOUT",&NOUT,PETSC_NULL);
170:   for (iout=1; iout<=NOUT; iout++){
171:     TSSetDuration(ts,time_steps,iout*ftime_original/NOUT);
172:     TSSolve(ts,global,&ftime);
173:     TSSetInitialTimeStep(ts,ftime,dt);
174:   }
175:   /* Interpolate solution at tfinal */
176:   TSGetSolution(ts,&global);
177:   TSInterpolate(ts,ftime_original,global);

179:   PetscOptionsHasName(PETSC_NULL,"-matlab_view",&flg);
180:   if (flg){ /* print solution into a MATLAB file */
181:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
182:     PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
183:     VecView(global,viewfile);
184:     PetscViewerDestroy(&viewfile);
185:   }

187:   /* display solver info for Sundials */
188:   TSGetType(ts,&tstype);
189:   PetscTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);
190:   if (sundials){
191:     PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
192:     TSView(ts,viewer);
193:     PetscViewerDestroy(&viewer);
194:     PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
195:     PCView(pc,viewer);
196:     PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",
197:                      size,tsinfo,pcinfo);
198:     PetscViewerDestroy(&viewer);
199:   }

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

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

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

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

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

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

242:   /* Initialize the array */
243:   VecGetArray(global,&localptr);

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

253:   VecRestoreArray(global,&localptr);
254:   return(0);
255: }

259: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
260: {
261:   VecScatter     scatter;
262:   IS             from,to;
263:   PetscInt       i,n,*idx,nsteps,maxsteps;
264:   Vec            tmp_vec;
266:   PetscScalar    *tmp;
267:   PetscReal      maxtime;
268:   Data           *data = (Data*)ctx;
269:   PetscReal      tfinal = data->tfinal;

272:   if (time > tfinal) return(0);

274:   TSGetTimeStepNumber(ts,&nsteps);
275:   /* display output at selected time steps */
276:   TSGetDuration(ts, &maxsteps, &maxtime);
277:   if (nsteps % 10 != 0 && time < maxtime) return(0);

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

282:   /* Set the index sets */
283:   PetscMalloc(n*sizeof(PetscInt),&idx);
284:   for(i=0; i<n; i++) idx[i]=i;

286:   /* Create local sequential vectors */
287:   VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);

289:   /* Create scatter context */
290:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
291:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
292:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
293:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
294:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

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

300:   PetscFree(idx);
301:   ISDestroy(&from);
302:   ISDestroy(&to);
303:   VecScatterDestroy(&scatter);
304:   VecDestroy(&tmp_vec);
305:   return(0);
306: }

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

321:   m = data->m;
322:   n = data->n;
323:   mn = m*n;
324:   dx = data->dx;
325:   dy = data->dy;
326:   a = data->a;
327:   epsilon = data->epsilon;

329:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
330:   xl = 0.5*a/dx+epsilon/(dx*dx);
331:   xr = -0.5*a/dx+epsilon/(dx*dx);
332:   yl = 0.5*a/dy+epsilon/(dy*dy);
333:   yr = -0.5*a/dy+epsilon/(dy*dy);

335:   row=0;
336:   v[0] = xc; v[1]=xr; v[2]=yr;
337:   idx[0]=0; idx[1]=2; idx[2]=m;
338:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

340:   row=m-1;
341:   v[0]=2.0*xl; v[1]=xc; v[2]=yr;
342:   idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
343:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

345:   for (i=1; i<m-1; i++) {
346:     row=i;
347:     v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
348:     idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
349:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
350:   }

352:   for (j=1; j<n-1; j++) {
353:     row=j*m;
354:     v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
355:     idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
356:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);

358:     row=j*m+m-1;
359:     v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
360:     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;
361:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);

363:     for(i=1; i<m-1; i++) {
364:       row=j*m+i;
365:       v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
366:       idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
367:       idx[4]=j*m+i+m;
368:       MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
369:     }
370:   }

372:   row=mn-m;
373:   v[0] = xc; v[1]=xr; v[2]=2.0*yl;
374:   idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
375:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

377:   row=mn-1;
378:   v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
379:   idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
380:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);

382:   for (i=1; i<m-1; i++) {
383:     row=mn-m+i;
384:     v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
385:     idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
386:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
387:   }

389:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
390:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

392:   /* *flag = SAME_NONZERO_PATTERN; */
393:   *flag = DIFFERENT_NONZERO_PATTERN;
394:   return(0);
395: }

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

416:   m       = data->m;
417:   n       = data->n;
418:   mn      = m*n;
419:   dx      = data->dx;
420:   dy      = data->dy;
421:   a       = data->a;
422:   epsilon = data->epsilon;

424:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
425:   xl = 0.5*a/dx+epsilon/(dx*dx);
426:   xr = -0.5*a/dx+epsilon/(dx*dx);
427:   yl = 0.5*a/dy+epsilon/(dy*dy);
428:   yr = -0.5*a/dy+epsilon/(dy*dy);

430:   /* Get the length of parallel vector */
431:   VecGetSize(globalin,&len);

433:   /* Set the index sets */
434:   PetscMalloc(len*sizeof(PetscInt),&idx);
435:   for(i=0; i<len; i++) idx[i]=i;

437:   /* Create local sequential vectors */
438:   VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
439:   VecDuplicate(tmp_in,&tmp_out);

441:   /* Create scatter context */
442:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
443:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
444:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
445:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
446:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
447:   VecScatterDestroy(&scatter);

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

452:   /* Extract outcome array*/
453:   VecGetArray(tmp_out,&outptr);

455:   outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
456:   outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
457:   for (i=1; i<m-1; i++) {
458:     outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
459:       +yr*inptr[i+m];
460:   }

462:   for (j=1; j<n-1; j++) {
463:     outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
464:       yl*inptr[j*m-m]+yr*inptr[j*m+m];
465:     outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
466:       yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
467:     for(i=1; i<m-1; i++) {
468:       outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
469:         +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
470:     }
471:   }

473:   outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
474:   outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
475:   for (i=1; i<m-1; i++) {
476:     outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
477:       +2*yl*inptr[mn-m+i-m];
478:   }

480:   VecRestoreArray(tmp_in,&inptr);
481:   VecRestoreArray(tmp_out,&outptr);

483:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
484:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
485:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);

487:   /* Destroy idx aand scatter */
488:   VecDestroy(&tmp_in);
489:   VecDestroy(&tmp_out);
490:   ISDestroy(&from);
491:   ISDestroy(&to);
492:   VecScatterDestroy(&scatter);

494:   PetscFree(idx);
495:   return(0);
496: }

500: PetscErrorCode PostStep(TS ts)
501: {
502:   PetscErrorCode  ierr;
503:   PetscReal       t;

506:   TSGetTime(ts,&t);
507:   PetscPrintf(PETSC_COMM_SELF,"  PostStep, t: %g\n",t);
508:   return(0);
509: }