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 petscsys.h
 19:  #include petscts.h


 24: typedef struct
 25: {
 26:   PetscInt         m;        /* the number of mesh points in x-direction */
 27:   PetscInt         n;      /* the number of mesh points in y-direction */
 28:   PetscReal         dx;     /* the grid space in x-direction */
 29:   PetscReal     dy;     /* the grid space in y-direction */
 30:   PetscReal     a;      /* the convection coefficient    */
 31:   PetscReal     epsilon; /* the diffusion coefficient    */
 32: } Data;

 34: /* two temporal functions */

 40: /* #undef PETSC_HAVE_SUNDIALS */

 44: int main(int argc,char **argv)
 45: {
 47:   PetscInt       time_steps = 100,steps;
 48:   PetscMPIInt    size;
 49:   Vec            global;
 50:   PetscReal      dt,ftime;
 51:   TS             ts;
 52:   PetscViewer         viewfile;
 53:   MatStructure   J_structure;
 54:   Mat            J = 0;
 55:   //SNES           snes;
 56:   Vec                  x;
 57:   Data                 data;
 58:   PetscInt          mn;
 59:   PetscTruth     flg;
 60: #if defined(PETSC_HAVE_SUNDIALS)
 61:   PC                 pc;
 62:   PetscViewer    viewer;
 63:   char           pcinfo[120],tsinfo[120];
 64: #endif

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

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

 91:   /* set user provided RHSFunction and RHSJacobian */
 92:   TSSetRHSFunction(ts,RHSFunction,&data);
 93:   MatCreate(PETSC_COMM_WORLD,&J);
 94:   MatSetSizes(J,PETSC_DECIDE,PETSC_DECIDE,mn,mn);
 95:   MatSetFromOptions(J);
 96:   RHSJacobian(ts,0.0,global,&J,&J,&J_structure,&data);
 97:   TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);

 99:   /* Use SUNDIALS */
100: #if defined(PETSC_HAVE_SUNDIALS)
101:   TSSetType(ts,TS_SUNDIALS);
102: #else
103:   TSSetType(ts,TS_EULER);
104: #endif
105:   dt   = 0.1;
106:   TSSetInitialTimeStep(ts,0.0,dt);
107:   TSSetDuration(ts,time_steps,1);
108:   TSSetSolution(ts,global);


111:   /* Pick up a Petsc preconditioner */
112:   /* one can always set method or preconditioner during the run time */
113: #if defined(PETSC_HAVE_SUNDIALS)
114:   TSSundialsGetPC(ts,&pc);
115:   PCSetType(pc,PCJACOBI);
116:   //TSSundialsSetType(ts,SUNDIALS_ADAMS);
117: #endif
118:   TSSetFromOptions(ts);

120:   TSSetUp(ts);
121:   TSStep(ts,&steps,&ftime);

123:   PetscOptionsHasName(PETSC_NULL,"-matlab_view",&flg);
124:   if (flg){ /* print solution into a MATLAB file */
125:     TSGetSolution(ts,&global);
126:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
127:     PetscViewerSetFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
128:     VecView(global,viewfile);
129:     PetscViewerDestroy(viewfile);
130:   }

132: #if defined(PETSC_HAVE_SUNDIALS)
133:   /* extracts the PC  from ts */
134:   TSSundialsGetPC(ts,&pc);
135:   PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
136:   TSView(ts,viewer);
137:   PetscViewerDestroy(viewer);
138:   PetscViewerStringOpen(PETSC_COMM_WORLD,pcinfo,120,&viewer);
139:   PCView(pc,viewer);
140:   PetscPrintf(PETSC_COMM_WORLD,"%d Procs,%s TSType, %s Preconditioner\n",
141:                      size,tsinfo,pcinfo);
142:   PetscViewerDestroy(viewer);
143: #endif

145:   /* free the memories */
146:   TSDestroy(ts);
147:   VecDestroy(global);
148:   VecDestroy(x);
149:   if (J) {ierr= MatDestroy(J);}
150:   //SNESDestroy(snes);
151:   PetscFinalize();
152:   return 0;
153: }

155: /* -------------------------------------------------------------------*/
156: /* the initial function */
157: PetscReal f_ini(PetscReal x,PetscReal y)
158: {
159:   PetscReal f;
160:   f=exp(-20.0*(pow(x-0.5,2.0)+pow(y-0.5,2.0)));
161:   return f;
162: }

166: PetscErrorCode Initial(Vec global,void *ctx)
167: {
168:   Data           *data = (Data*)ctx;
169:   PetscInt       m,row,col;
170:   PetscReal      x,y,dx,dy;
171:   PetscScalar    *localptr;
172:   PetscInt       i,mybase,myend,locsize;

175:   /* make the local  copies of parameters */
176:   m = data->m;
177:   dx = data->dx;
178:   dy = data->dy;

180:   /* determine starting point of each processor */
181:   VecGetOwnershipRange(global,&mybase,&myend);
182:   VecGetLocalSize(global,&locsize);

184:   /* Initialize the array */
185:   VecGetArray(global,&localptr);

187:   for (i=0; i<locsize; i++) {
188:     row = 1+(mybase+i)-((mybase+i)/m)*m;
189:     col = (mybase+i)/m+1;
190:     x = dx*row;
191:     y = dy*col;
192:     localptr[i] = f_ini(x,y);
193:   }
194: 
195:   VecRestoreArray(global,&localptr);
196:   return 0;
197: }

201: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
202: {
203:   VecScatter     scatter;
204:   IS             from,to;
205:   PetscInt       i,n,*idx;
206:   Vec            tmp_vec;
208:   PetscScalar    *tmp;

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

213:   /* Set the index sets */
214:   PetscMalloc(n*sizeof(PetscInt),&idx);
215:   for(i=0; i<n; i++) idx[i]=i;
216: 
217:   /* Create local sequential vectors */
218:   VecCreateSeq(PETSC_COMM_SELF,n,&tmp_vec);

220:   /* Create scatter context */
221:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&from);
222:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,&to);
223:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
224:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
225:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

227:   VecGetArray(tmp_vec,&tmp);
228:   PetscPrintf(PETSC_COMM_WORLD,"At t =%14.6e u= %14.6e at the center \n",time,PetscRealPart(tmp[n/2]));
229:   VecRestoreArray(tmp_vec,&tmp);

231:   PetscFree(idx);
232:   ISDestroy(from);
233:   ISDestroy(to);
234:   VecScatterDestroy(scatter);
235:   VecDestroy(tmp_vec);
236:   return 0;
237: }

239: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
242: PetscErrorCode FormFunction(SNES snes,Vec globalin,Vec globalout,void *ptr)
243: {
244:   Data           *data = (Data*)ptr;
245:   PetscInt       m,n,mn;
246:   PetscReal      dx,dy;
247:   PetscReal      xc,xl,xr,yl,yr;
248:   PetscReal      a,epsilon;
249:   PetscScalar    *inptr,*outptr;
250:   PetscInt       i,j,len;
252:   IS             from,to;
253:   PetscInt       *idx;
254:   VecScatter     scatter;
255:   Vec            tmp_in,tmp_out;

257:   m       = data->m;
258:   n       = data->n;
259:   mn      = m*n;
260:   dx      = data->dx;
261:   dy      = data->dy;
262:   a       = data->a;
263:   epsilon = data->epsilon;

265:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
266:   xl = 0.5*a/dx+epsilon/(dx*dx);
267:   xr = -0.5*a/dx+epsilon/(dx*dx);
268:   yl = 0.5*a/dy+epsilon/(dy*dy);
269:   yr = -0.5*a/dy+epsilon/(dy*dy);
270: 
271:   /* Get the length of parallel vector */
272:   VecGetSize(globalin,&len);

274:   /* Set the index sets */
275:   PetscMalloc(len*sizeof(PetscInt),&idx);
276:   for(i=0; i<len; i++) idx[i]=i;
277: 
278:   /* Create local sequential vectors */
279:   VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
280:   VecDuplicate(tmp_in,&tmp_out);

282:   /* Create scatter context */
283:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,&from);
284:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,&to);
285:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
286:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
287:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
288:   VecScatterDestroy(scatter);

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

293:   /* Extract outcome array*/
294:   VecGetArray(tmp_out,&outptr);

296:   outptr[0] = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
297:   outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
298:   for (i=1; i<m-1; i++) {
299:     outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]
300:       +yr*inptr[i+m];
301:   }

303:   for (j=1; j<n-1; j++) {
304:     outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+
305:       yl*inptr[j*m-m]+yr*inptr[j*m+m];
306:     outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+
307:       yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
308:     for(i=1; i<m-1; i++) {
309:       outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]
310:         +yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
311:     }
312:   }

314:   outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
315:   outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
316:   for (i=1; i<m-1; i++) {
317:     outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]
318:       +2*yl*inptr[mn-m+i-m];
319:   }

321:   VecRestoreArray(tmp_in,&inptr);
322:   VecRestoreArray(tmp_out,&outptr);

324:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
325:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
326:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
327: 
328:   /* Destroy idx aand scatter */
329:   VecDestroy(tmp_in);
330:   VecDestroy(tmp_out);
331:   ISDestroy(from);
332:   ISDestroy(to);
333:   VecScatterDestroy(scatter);

335:   PetscFree(idx);
336:   return 0;
337: }

341: PetscErrorCode FormJacobian(SNES snes,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
342: {
343:   Data           *data = (Data*)ptr;
344:   Mat            A = *AA;
345:   PetscScalar    v[1],one = 1.0;
346:   PetscInt       idx[1],i,j,row;
348:   PetscInt       m,n,mn;

350:   m = data->m;
351:   n = data->n;
352:   mn = (data->m)*(data->n);
353: 
354:   for(i=0; i<mn; i++) {
355:     idx[0] = i;
356:     v[0] = one;
357:     MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
358:   }

360:   for(i=0; i<mn-m; i++) {
361:     idx[0] = i+m;
362:     v[0]= one;
363:     MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
364:   }

366:   for(i=m; i<mn; i++) {
367:     idx[0] = i-m;
368:     v[0] = one;
369:     MatSetValues(A,1,&i,1,idx,v,INSERT_VALUES);
370:   }

372:   for(j=0; j<n; j++) {
373:     for(i=0; i<m-1; i++) {
374:       row = j*m+i;
375:       idx[0] = j*m+i+1;
376:       v[0] = one;
377:       MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
378:     }
379:   }

381:   for(j=0; j<n; j++) {
382:     for(i=1; i<m; i++) {
383:       row = j*m+i;
384:       idx[0] = j*m+i-1;
385:       v[0] = one;
386:       MatSetValues(A,1,&row,1,idx,v,INSERT_VALUES);
387:     }
388:   }
389:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
390:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
391:   *flag = SAME_NONZERO_PATTERN;
392:   return 0;
393: }

397: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *flag,void *ptr)
398: {
399:   Data           *data = (Data*)ptr;
400:   Mat            A = *AA;
401:   PetscScalar    v[5];
402:   PetscInt       idx[5],i,j,row;
404:   PetscInt       m,n,mn;
405:   PetscReal      dx,dy,a,epsilon,xc,xl,xr,yl,yr;

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

415:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
416:   xl = 0.5*a/dx+epsilon/(dx*dx);
417:   xr = -0.5*a/dx+epsilon/(dx*dx);
418:   yl = 0.5*a/dy+epsilon/(dy*dy);
419:   yr = -0.5*a/dy+epsilon/(dy*dy);

421:   row=0;
422:   v[0] = xc; v[1]=xr; v[2]=yr;
423:   idx[0]=0; idx[1]=2; idx[2]=m;
424:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

426:   row=m-1;
427:   v[0]=2.0*xl; v[1]=xc; v[2]=yr;
428:   idx[0]=m-2; idx[1]=m-1; idx[2]=m-1+m;
429:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

431:   for (i=1; i<m-1; i++) {
432:     row=i;
433:     v[0]=xl; v[1]=xc; v[2]=xr; v[3]=yr;
434:     idx[0]=i-1; idx[1]=i; idx[2]=i+1; idx[3]=i+m;
435:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
436:   }

438:   for (j=1; j<n-1; j++) {
439:     row=j*m;
440:     v[0]=xc; v[1]=xr; v[2]=yl; v[3]=yr;
441:     idx[0]=j*m; idx[1]=j*m; idx[2]=j*m-m; idx[3]=j*m+m;
442:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
443: 
444:     row=j*m+m-1;
445:     v[0]=xc; v[1]=2.0*xl; v[2]=yl; v[3]=yr;
446:     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;
447:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);

449:     for(i=1; i<m-1; i++) {
450:       row=j*m+i;
451:       v[0]=xc; v[1]=xl; v[2]=xr; v[3]=yl; v[4]=yr;
452:       idx[0]=j*m+i; idx[1]=j*m+i-1; idx[2]=j*m+i+1; idx[3]=j*m+i-m;
453:       idx[4]=j*m+i+m;
454:       MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
455:     }
456:   }

458:   row=mn-m;
459:   v[0] = xc; v[1]=xr; v[2]=2.0*yl;
460:   idx[0]=mn-m; idx[1]=mn-m+1; idx[2]=mn-m-m;
461:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);
462: 
463:   row=mn-1;
464:   v[0] = xc; v[1]=2.0*xl; v[2]=2.0*yl;
465:   idx[0]=mn-1; idx[1]=mn-2; idx[2]=mn-1-m;
466:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);

468:   for (i=1; i<m-1; i++) {
469:     row=mn-m+i;
470:     v[0]=xl; v[1]=xc; v[2]=xr; v[3]=2.0*yl;
471:     idx[0]=mn-m+i-1; idx[1]=mn-m+i; idx[2]=mn-m+i+1; idx[3]=mn-m+i-m;
472:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
473:   }


476:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
477:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

479:   /* *flag = SAME_NONZERO_PATTERN; */
480:   *flag = DIFFERENT_NONZERO_PATTERN;
481:   return 0;
482: }

486: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
487: {
489:   SNES           snes = PETSC_NULL;

491:   FormFunction(snes,globalin,globalout,ctx);
492:   return 0;
493: }