Actual source code: ex4.c

petsc-3.9.4 2018-09-11
Report Typos and Errors
  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: {
 40:   PetscInt       time_steps=100,iout,NOUT=1;
 41:   PetscMPIInt    size;
 42:   Vec            global;
 43:   PetscReal      dt,ftime,ftime_original;
 44:   TS             ts;
 45:   PetscViewer    viewfile;
 46:   Mat            J = 0;
 47:   Vec            x;
 48:   Data           data;
 49:   PetscInt       mn;
 50:   PetscBool      flg;
 51:   MatColoring    mc;
 52:   ISColoring     iscoloring;
 53:   MatFDColoring  matfdcoloring        = 0;
 54:   PetscBool      fd_jacobian_coloring = PETSC_FALSE;
 55:   SNES           snes;
 56:   KSP            ksp;
 57:   PC             pc;

 59:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
 60:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

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

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

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

 86:   TSSetTimeStep(ts,dt);
 87:   TSSetMaxSteps(ts,time_steps);
 88:   TSSetMaxTime(ts,ftime_original);
 89:   TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
 90:   TSSetSolution(ts,global);

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

100:   PetscOptionsHasName(NULL,NULL,"-ts_fd",&flg);
101:   if (!flg) {
102:     TSSetRHSJacobian(ts,J,J,RHSJacobian,&data);
103:   } else {
104:     TSGetSNES(ts,&snes);
105:     PetscOptionsHasName(NULL,NULL,"-fd_color",&fd_jacobian_coloring);
106:     if (fd_jacobian_coloring) { /* Use finite differences with coloring */
107:       /* Get data structure of J */
108:       PetscBool pc_diagonal;
109:       PetscOptionsHasName(NULL,NULL,"-pc_diagonal",&pc_diagonal);
110:       if (pc_diagonal) { /* the preconditioner of J is a diagonal matrix */
111:         PetscInt    rstart,rend,i;
112:         PetscScalar zero=0.0;
113:         MatGetOwnershipRange(J,&rstart,&rend);
114:         for (i=rstart; i<rend; i++) {
115:           MatSetValues(J,1,&i,1,&i,&zero,INSERT_VALUES);
116:         }
117:         MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
118:         MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
119:       } else {
120:         /* Fill the structure using the expensive SNESComputeJacobianDefault. Temporarily set up the TS so we can call this function */
121:         TSSetType(ts,TSBEULER);
122:         TSSetUp(ts);
123:         SNESComputeJacobianDefault(snes,x,J,J,ts);
124:       }

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

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

151:   TSSetFromOptions(ts);
152:   TSSetUp(ts);

154:   /* Test TSSetPostStep() */
155:   PetscOptionsHasName(NULL,NULL,"-test_PostStep",&flg);
156:   if (flg) {
157:     TSSetPostStep(ts,PostStep);
158:   }

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

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

182:   /* free the memories */
183:   TSDestroy(&ts);
184:   VecDestroy(&global);
185:   VecDestroy(&x);
186:   MatDestroy(&J);
187:   if (fd_jacobian_coloring) {MatFDColoringDestroy(&matfdcoloring);}
188:   PetscFinalize();
189:   return ierr;
190: }

192: /* -------------------------------------------------------------------*/
193: /* the initial function */
194: PetscReal f_ini(PetscReal x,PetscReal y)
195: {
196:   PetscReal f;

198:   f=PetscExpReal(-20.0*(PetscPowRealInt(x-0.5,2)+PetscPowRealInt(y-0.5,2)));
199:   return f;
200: }

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

212:   /* make the local  copies of parameters */
213:   m  = data->m;
214:   dx = data->dx;
215:   dy = data->dy;

217:   /* determine starting point of each processor */
218:   VecGetOwnershipRange(global,&mybase,&myend);
219:   VecGetLocalSize(global,&locsize);

221:   /* Initialize the array */
222:   VecGetArray(global,&localptr);

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

232:   VecRestoreArray(global,&localptr);
233:   return(0);
234: }

236: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec global,void *ctx)
237: {
238:   VecScatter        scatter;
239:   IS                from,to;
240:   PetscInt          i,n,*idx,nsteps,maxsteps;
241:   Vec               tmp_vec;
242:   PetscErrorCode    ierr;
243:   const PetscScalar *tmp;

246:   TSGetStepNumber(ts,&nsteps);
247:   /* display output at selected time steps */
248:   TSGetMaxSteps(ts, &maxsteps);
249:   if (nsteps % 10 != 0) return(0);

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

254:   /* Set the index sets */
255:   PetscMalloc1(n,&idx);
256:   for (i=0; i<n; i++) idx[i]=i;

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

261:   /* Create scatter context */
262:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&from);
263:   ISCreateGeneral(PETSC_COMM_SELF,n,idx,PETSC_COPY_VALUES,&to);
264:   VecScatterCreate(global,from,tmp_vec,to,&scatter);
265:   VecScatterBegin(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);
266:   VecScatterEnd(scatter,global,tmp_vec,INSERT_VALUES,SCATTER_FORWARD);

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

272:   PetscFree(idx);
273:   ISDestroy(&from);
274:   ISDestroy(&to);
275:   VecScatterDestroy(&scatter);
276:   VecDestroy(&tmp_vec);
277:   return(0);
278: }

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

290:   m       = data->m;
291:   n       = data->n;
292:   mn      = m*n;
293:   dx      = data->dx;
294:   dy      = data->dy;
295:   a       = data->a;
296:   epsilon = data->epsilon;

298:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
299:   xl = 0.5*a/dx+epsilon/(dx*dx);
300:   xr = -0.5*a/dx+epsilon/(dx*dx);
301:   yl = 0.5*a/dy+epsilon/(dy*dy);
302:   yr = -0.5*a/dy+epsilon/(dy*dy);

304:   row    = 0;
305:   v[0]   = xc;  v[1] = xr;  v[2] = yr;
306:   idx[0] = 0; idx[1] = 2; idx[2] = m;
307:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

309:   row    = m-1;
310:   v[0]   = 2.0*xl; v[1] = xc;    v[2] = yr;
311:   idx[0] = m-2;  idx[1] = m-1; idx[2] = m-1+m;
312:   MatSetValues(A,1,&row,3,idx,v,INSERT_VALUES);

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

321:   for (j=1; j<n-1; j++) {
322:     row    = j*m;
323:     v[0]   = xc;    v[1] = xr;    v[2] = yl;      v[3] = yr;
324:     idx[0] = j*m; idx[1] = j*m; idx[2] = j*m-m; idx[3] = j*m+m;
325:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);

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

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

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

346:   row    = mn-1;
347:   v[0]   = xc;     v[1] = 2.0*xl; v[2] = 2.0*yl;
348:   idx[0] = mn-1; idx[1] = mn-2; idx[2] = mn-1-m;
349:   MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);

351:   for (i=1; i<m-1; i++) {
352:     row    = mn-m+i;
353:     v[0]   = xl;         v[1] = xc;       v[2] = xr;         v[3] = 2.0*yl;
354:     idx[0] = mn-m+i-1; idx[1] = mn-m+i; idx[2] = mn-m+i+1; idx[3] = mn-m+i-m;
355:     MatSetValues(A,1,&row,4,idx,v,INSERT_VALUES);
356:   }

358:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
359:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);

361:   return(0);
362: }

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

382:   m       = data->m;
383:   n       = data->n;
384:   mn      = m*n;
385:   dx      = data->dx;
386:   dy      = data->dy;
387:   a       = data->a;
388:   epsilon = data->epsilon;

390:   xc = -2.0*epsilon*(1.0/(dx*dx)+1.0/(dy*dy));
391:   xl = 0.5*a/dx+epsilon/(dx*dx);
392:   xr = -0.5*a/dx+epsilon/(dx*dx);
393:   yl = 0.5*a/dy+epsilon/(dy*dy);
394:   yr = -0.5*a/dy+epsilon/(dy*dy);

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

399:   /* Set the index sets */
400:   PetscMalloc1(len,&idx);
401:   for (i=0; i<len; i++) idx[i]=i;

403:   /* Create local sequential vectors */
404:   VecCreateSeq(PETSC_COMM_SELF,len,&tmp_in);
405:   VecDuplicate(tmp_in,&tmp_out);

407:   /* Create scatter context */
408:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&from);
409:   ISCreateGeneral(PETSC_COMM_SELF,len,idx,PETSC_COPY_VALUES,&to);
410:   VecScatterCreate(globalin,from,tmp_in,to,&scatter);
411:   VecScatterBegin(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
412:   VecScatterEnd(scatter,globalin,tmp_in,INSERT_VALUES,SCATTER_FORWARD);
413:   VecScatterDestroy(&scatter);

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

418:   /* Extract outcome array*/
419:   VecGetArray(tmp_out,&outptr);

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

427:   for (j=1; j<n-1; j++) {
428:     outptr[j*m] = xc*inptr[j*m]+xr*inptr[j*m+1]+ yl*inptr[j*m-m]+yr*inptr[j*m+m];
429:     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];
430:     for (i=1; i<m-1; i++) {
431:       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];
432:     }
433:   }

435:   outptr[mn-m] = xc*inptr[mn-m]+xr*inptr[mn-m+1]+2.0*yl*inptr[mn-m-m];
436:   outptr[mn-1] = 2.0*xl*inptr[mn-2]+xc*inptr[mn-1]+2.0*yl*inptr[mn-1-m];
437:   for (i=1; i<m-1; i++) {
438:     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];
439:   }

441:   VecRestoreArrayRead(tmp_in,&inptr);
442:   VecRestoreArray(tmp_out,&outptr);

444:   VecScatterCreate(tmp_out,from,globalout,to,&scatter);
445:   VecScatterBegin(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);
446:   VecScatterEnd(scatter,tmp_out,globalout,INSERT_VALUES,SCATTER_FORWARD);

448:   /* Destroy idx aand scatter */
449:   VecDestroy(&tmp_in);
450:   VecDestroy(&tmp_out);
451:   ISDestroy(&from);
452:   ISDestroy(&to);
453:   VecScatterDestroy(&scatter);

455:   PetscFree(idx);
456:   return(0);
457: }

459: PetscErrorCode PostStep(TS ts)
460: {
462:   PetscReal      t;

465:   TSGetTime(ts,&t);
466:   PetscPrintf(PETSC_COMM_SELF,"  PostStep, t: %g\n",(double)t);
467:   return(0);
468: }

470: /*TEST

472:     test:
473:       args: -ts_fd -ts_type beuler
474:       output_file: output/ex4.out

476:     test:
477:       suffix: 2
478:       args: -ts_fd -ts_type beuler
479:       nsize: 2
480:       output_file: output/ex4.out

482:     test:
483:       suffix: 3
484:       args: -ts_fd -ts_type cn

486:     test:
487:       suffix: 4
488:       args: -ts_fd -ts_type cn
489:       output_file: output/ex4_3.out
490:       nsize: 2

492:     test:
493:       suffix: 5
494:       args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
495:       output_file: output/ex4.out

497:     test:
498:       suffix: 6
499:       args: -ts_type beuler -ts_fd -fd_color -mat_coloring_type sl
500:       output_file: output/ex4.out
501:       nsize: 2

503:     test:
504:       suffix: 7
505:       requires: !single
506:       args: -ts_fd -ts_type beuler -test_PostStep -ts_dt .1

508: TEST*/