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.

 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,void*);
 35: extern PetscErrorCode PostStep(TS);

 37: int main(int argc,char **argv)
 38: {
 39:   PetscInt       time_steps=100,iout,NOUT=1;
 40:   Vec            global;
 41:   PetscReal      dt,ftime,ftime_original;
 42:   TS             ts;
 43:   PetscViewer    viewfile;
 44:   Mat            J = 0;
 45:   Vec            x;
 46:   Data           data;
 47:   PetscInt       mn;
 48:   PetscBool      flg;
 49:   MatColoring    mc;
 50:   ISColoring     iscoloring;
 51:   MatFDColoring  matfdcoloring        = 0;
 52:   PetscBool      fd_jacobian_coloring = PETSC_FALSE;
 53:   SNES           snes;
 54:   KSP            ksp;
 55:   PC             pc;

 57:   PetscInitialize(&argc,&argv,(char*)0,help);

 59:   /* set data */
 60:   data.m       = 9;
 61:   data.n       = 9;
 62:   data.a       = 1.0;
 63:   data.epsilon = 0.1;
 64:   data.dx      = 1.0/(data.m+1.0);
 65:   data.dy      = 1.0/(data.n+1.0);
 66:   mn           = (data.m)*(data.n);
 67:   PetscOptionsGetInt(NULL,NULL,"-time",&time_steps,NULL);

 69:   /* set initial conditions */
 70:   VecCreate(PETSC_COMM_WORLD,&global);
 71:   VecSetSizes(global,PETSC_DECIDE,mn);
 72:   VecSetFromOptions(global);
 73:   Initial(global,&data);
 74:   VecDuplicate(global,&x);

 76:   /* create timestep context */
 77:   TSCreate(PETSC_COMM_WORLD,&ts);
 78:   TSMonitorSet(ts,Monitor,&data,NULL);
 79:   TSSetType(ts,TSEULER);
 80:   dt   = 0.1;
 81:   ftime_original = data.tfinal = 1.0;

 83:   TSSetTimeStep(ts,dt);
 84:   TSSetMaxSteps(ts,time_steps);
 85:   TSSetMaxTime(ts,ftime_original);
 86:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
 87:   TSSetSolution(ts,global);

 89:   /* set user provided RHSFunction and RHSJacobian */
 90:   TSSetRHSFunction(ts,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,NULL);
 95:   MatMPIAIJSetPreallocation(J,5,NULL,5,NULL);

 97:   PetscOptionsHasName(NULL,NULL,"-ts_fd",&flg);
 98:   if (!flg) {
 99:     TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
100:   } else {
101:     TSGetSNES(ts,&snes);
102:     PetscOptionsHasName(NULL,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:       PetscOptionsHasName(NULL,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 {
117:         /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
118:         TSSetType(ts,TSBEULER);
119:         TSSetUp(ts);
120:         SNESComputeJacobianDefault(snes,x,J,J,ts);
121:       }

123:       /* create coloring context */
124:       MatColoringCreate(J,&mc);
125:       MatColoringSetType(mc,MATCOLORINGSL);
126:       MatColoringSetFromOptions(mc);
127:       MatColoringApply(mc,&iscoloring);
128:       MatColoringDestroy(&mc);
129:       MatFDColoringCreate(J,iscoloring,&matfdcoloring);
130:       MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
131:       MatFDColoringSetFromOptions(matfdcoloring);
132:       MatFDColoringSetUp(J,iscoloring,matfdcoloring);
133:       SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);
134:       ISColoringDestroy(&iscoloring);
135:     } else { /* Use finite differences (slow) */
136:       SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);
137:     }
138:   }

140:   /* Pick up a Petsc preconditioner */
141:   /* one can always set method or preconditioner during the run time */
142:   TSGetSNES(ts,&snes);
143:   SNESGetKSP(snes,&ksp);
144:   KSPGetPC(ksp,&pc);
145:   PCSetType(pc,PCJACOBI);
146:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);

148:   TSSetFromOptions(ts);
149:   TSSetUp(ts);

151:   /* Test TSSetPostStep() */
152:   PetscOptionsHasName(NULL,NULL,"-test_PostStep",&flg);
153:   if (flg) {
154:     TSSetPostStep(ts,PostStep);
155:   }

157:   PetscOptionsGetInt(NULL,NULL,"-NOUT",&NOUT,NULL);
158:   for (iout=1; iout<=NOUT; iout++) {
159:     TSSetMaxSteps(ts,time_steps);
160:     TSSetMaxTime(ts,iout*ftime_original/NOUT);
161:     TSSolve(ts,global);
162:     TSGetSolveTime(ts,&ftime);
163:     TSSetTime(ts,ftime);
164:     TSSetTimeStep(ts,dt);
165:   }
166:   /* Interpolate solution at tfinal */
167:   TSGetSolution(ts,&global);
168:   TSInterpolate(ts,ftime_original,global);

170:   PetscOptionsHasName(NULL,NULL,"-matlab_view",&flg);
171:   if (flg) { /* print solution into a MATLAB file */
172:     PetscViewerASCIIOpen(PETSC_COMM_WORLD,"out.m",&viewfile);
173:     PetscViewerPushFormat(viewfile,PETSC_VIEWER_ASCII_MATLAB);
174:     VecView(global,viewfile);
175:     PetscViewerPopFormat(viewfile);
176:     PetscViewerDestroy(&viewfile);
177:   }

179:   /* free the memories */
180:   TSDestroy(&ts);
181:   VecDestroy(&global);
182:   VecDestroy(&x);
183:   MatDestroy(&J);
184:   if (fd_jacobian_coloring) MatFDColoringDestroy(&matfdcoloring);
185:   PetscFinalize();
186:   return 0;
187: }

189: /* -------------------------------------------------------------------*/
190: /* the initial function */
191: PetscReal f_ini(PetscReal x,PetscReal y)
192: {
193:   PetscReal f;

195:   f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2)));
196:   return f;
197: }

199: PetscErrorCode Initial(Vec global,void *ctx)
200: {
201:   Data           *data = (Data*)ctx;
202:   PetscInt       m,row,col;
203:   PetscReal      x,y,dx,dy;
204:   PetscScalar    *localptr;
205:   PetscInt       i,mybase,myend,locsize;

208:   /* make the local  copies of parameters */
209:   m  = data->m;
210:   dx = data->dx;
211:   dy = data->dy;

213:   /* determine starting point of each processor */
214:   VecGetOwnershipRange(global,&mybase,&myend);
215:   VecGetLocalSize(global,&locsize);

217:   /* Initialize the array */
218:   VecGetArrayWrite(global,&localptr);

220:   for (i=0; i<locsize; i++) {
221:     row         = 1+(mybase+i)-((mybase+i)/m)*m;
222:     col         = (mybase+i)/m+1;
223:     x           = dx*row;
224:     y           = dy*col;
225:     localptr[i] = f_ini(x,y);
226:   }

228:   VecRestoreArrayWrite(global,&localptr);
229:   return 0;
230: }

232: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
233: {
234:   VecScatter        scatter;
235:   IS                from,to;
236:   PetscInt          i,n,*idx,nsteps,maxsteps;
237:   Vec               tmp_vec;
238:   const PetscScalar *tmp;

241:   TSGetStepNumber(ts,&nsteps);
242:   /* display output at selected time steps */
243:   TSGetMaxSteps(ts, &maxsteps);
244:   if (nsteps % 10 != 0) return 0;

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

249:   /* Set the index sets */
250:   PetscMalloc1(n,&idx);
251:   for (i=0; i<n; i++) idx[i]=i;

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

256:   /* Create scatter context */
257:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
258:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
259:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
260:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
261:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

263:   VecGetArrayRead(tmp_vec,&tmp);
264:   PetscPrintf(PETSC_COMM_WORLD,"At t[%D] =%14.2e u= %14.2e at the center \n",nsteps,(double)time,(double)PetscRealPart(tmp[n/2]));
265:   VecRestoreArrayRead(tmp_vec,&tmp);

267:   PetscFree(idx);
268:   ISDestroy(&from);
269:   ISDestroy(&to);
270:   VecScatterDestroy(&scatter);
271:   VecDestroy(&tmp_vec);
272:   return 0;
273: }

275: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat BB,void *ptr)
276: {
277:   Data           *data = (Data*)ptr;
278:   PetscScalar    v[5];
279:   PetscInt       idx[5],i,j,row;
280:   PetscInt       m,n,mn;
281:   PetscReal      dx,dy,a,epsilon,xc,xl,xr,yl,yr;

284:   m       = data->m;
285:   n       = data->n;
286:   mn      = m*n;
287:   dx      = data->dx;
288:   dy      = data->dy;
289:   a       = data->a;
290:   epsilon = data->epsilon;

292:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
293:   xl = 0.5*a/dx+epsilon/(dx*dx);
294:   xr = -0.5*a/dx+epsilon/(dx*dx);
295:   yl = 0.5*a/dy+epsilon/(dy*dy);
296:   yr = -0.5*a/dy+epsilon/(dy*dy);

298:   row    = 0;
299:   v[0]   = xc;  v[1] = xr;  v[2] = yr;
300:   idx[0] = 0; idx[1] = 2; idx[2] = m;
301:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

303:   row    = m-1;
304:   v[0]   = 2.0*xl; v[1] = xc;    v[2] = yr;
305:   idx[0] = m-2;  idx[1] = m-1; idx[2] = m-1+m;
306:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

308:   for (i=1; i<m-1; i++) {
309:     row    = i;
310:     v[0]   = xl;    v[1] = xc;  v[2] = xr;    v[3] = yr;
311:     idx[0] = i-1; idx[1] = i; idx[2] = i+1; idx[3] = i+m;
312:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
313:   }

315:   for (j=1; j<n-1; j++) {
316:     row    = j*m;
317:     v[0]   = xc;    v[1] = xr;    v[2] = yl;      v[3] = yr;
318:     idx[0] = j*m; idx[1] = j*m; idx[2] = j*m-m; idx[3] = j*m+m;
319:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);

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

326:     for (i=1; i<m-1; i++) {
327:       row    = j*m+i;
328:       v[0]   = xc;      v[1] = xl;        v[2] = xr;        v[3] = yl; v[4]=yr;
329:       idx[0] = j*m+i; idx[1] = j*m+i-1; idx[2] = j*m+i+1; idx[3] = j*m+i-m;
330:       idx[4] = j*m+i+m;
331:       MatSetValues(A,1,&row,5,idx,v,INSERT_VALUES);
332:     }
333:   }

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

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

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

352:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
353:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

355:   return 0;
356: }

358: /* globalout = -a*(u_x+u_y) + epsilon*(u_xx+u_yy) */
359: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
360: {
361:   Data              *data = (Data*)ctx;
362:   PetscInt          m,n,mn;
363:   PetscReal         dx,dy;
364:   PetscReal         xc,xl,xr,yl,yr;
365:   PetscReal         a,epsilon;
366:   PetscScalar       *outptr;
367:   const PetscScalar *inptr;
368:   PetscInt          i,j,len;
369:   IS                from,to;
370:   PetscInt          *idx;
371:   VecScatter        scatter;
372:   Vec               tmp_in,tmp_out;

375:   m       = data->m;
376:   n       = data->n;
377:   mn      = m*n;
378:   dx      = data->dx;
379:   dy      = data->dy;
380:   a       = data->a;
381:   epsilon = data->epsilon;

383:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
384:   xl = 0.5*a/dx+epsilon/(dx*dx);
385:   xr = -0.5*a/dx+epsilon/(dx*dx);
386:   yl = 0.5*a/dy+epsilon/(dy*dy);
387:   yr = -0.5*a/dy+epsilon/(dy*dy);

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

392:   /* Set the index sets */
393:   PetscMalloc1(len,&idx);
394:   for (i=0; i<len; i++) idx[i]=i;

396:   /* Create local sequential vectors */
397:   VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
398:   VecDuplicate(tmp_in,&tmp_out);

400:   /* Create scatter context */
401:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
402:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
403:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
404:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
405:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
406:   VecScatterDestroy(&scatter);

408:   /*Extract income array - include ghost points */
409:   VecGetArrayRead(tmp_in,&inptr);

411:   /* Extract outcome array*/
412:   VecGetArrayWrite(tmp_out,&outptr);

414:   outptr[0]   = xc*inptr[0]+xr*inptr[1]+yr*inptr[m];
415:   outptr[m-1] = 2.0*xl*inptr[m-2]+xc*inptr[m-1]+yr*inptr[m-1+m];
416:   for (i=1; i<m-1; i++) {
417:     outptr[i] = xc*inptr[i]+xl*inptr[i-1]+xr*inptr[i+1]+yr*inptr[i+m];
418:   }

420:   for (j=1; j<n-1; j++) {
421:     outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+ yl*inptr[j*m-m]+yr*inptr[j*m+m];
422:     outptr[j*m+m-1] = xc*inptr[j*m+m-1]+2.0*xl*inptr[j*m+m-1-1]+ yl*inptr[j*m+m-1-m]+yr*inptr[j*m+m-1+m];
423:     for (i=1; i<m-1; i++) {
424:       outptr[j*m+i] = xc*inptr[j*m+i]+xl*inptr[j*m+i-1]+xr*inptr[j*m+i+1]+yl*inptr[j*m+i-m]+yr*inptr[j*m+i+m];
425:     }
426:   }

428:   outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
429:   outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
430:   for (i=1; i<m-1; i++) {
431:     outptr[mn-m+i] = xc*inptr[mn-m+i]+xl*inptr[mn-m+i-1]+xr*inptr[mn-m+i+1]+2*yl*inptr[mn-m+i-m];
432:   }

434:   VecRestoreArrayRead(tmp_in,&inptr);
435:   VecRestoreArrayWrite(tmp_out,&outptr);

437:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
438:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
439:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);

441:   /* Destroy idx aand scatter */
442:   VecDestroy(&tmp_in);
443:   VecDestroy(&tmp_out);
444:   ISDestroy(&from);
445:   ISDestroy(&to);
446:   VecScatterDestroy(&scatter);

448:   PetscFree(idx);
449:   return 0;
450: }

452: PetscErrorCode PostStep(TS ts)
453: {
454:   PetscReal      t;

457:   TSGetTime(ts,&t);
458:   PetscPrintf(PETSC_COMM_SELF,"  PostStep, t: %g\n",(double)t);
459:   return 0;
460: }

462: /*TEST

464:     test:
465:       args: -ts_fd -ts_type beuler
466:       output_file: output/ex4.out

468:     test:
469:       suffix: 2
470:       args: -ts_fd -ts_type beuler
471:       nsize: 2
472:       output_file: output/ex4.out

474:     test:
475:       suffix: 3
476:       args: -ts_fd -ts_type cn

478:     test:
479:       suffix: 4
480:       args: -ts_fd -ts_type cn
481:       output_file: output/ex4_3.out
482:       nsize: 2

484:     test:
485:       suffix: 5
486:       args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
487:       output_file: output/ex4.out

489:     test:
490:       suffix: 6
491:       args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
492:       output_file: output/ex4.out
493:       nsize: 2

495:     test:
496:       suffix: 7
497:       requires: !single
498:       args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1

500:     test:
501:       suffix: 8
502:       requires: !single
503:       args: -ts_type rk -ts_rk_type 5dp -ts_dt .01 -ts_adapt_type none -ts_view

505: TEST*/