Actual source code: ex1.c

  1: /*
  2:        Formatted test for TS routines.

  4:           Solves U_t = U_xx 
  5:      F(t,u) = (u_i+1 - 2u_i + u_i-1)/h^2
  6:        using several different schemes. 
  7: */

  9: /* Usage: 
 10:    ./ex1 -nox -ts_type beuler -ts_view 
 11:    ./ex1 -nox -linear_constant_matrix -ts_type beuler
 12:    ./ex1 -nox -linear_variable_matrix -ts_type beuler
 13: */

 15: static char help[] = "Solves 1D heat equation.\n\n";

 17:  #include petscda.h
 18:  #include petscsys.h
 19:  #include petscts.h

 21: #define PETSC_NEAR(a,b,c) (!(PetscAbsReal((a)-(b)) > (c)*PetscMax(PetscAbsReal(a),PetscAbsReal(b))))

 23: typedef struct {
 24:   Vec         global,local,localwork,solution;    /* location for local work (with ghost points) vector */
 25:   DA          da;                    /* manages ghost point communication */
 26:   PetscViewer viewer1,viewer2;
 27:   PetscInt    M;                     /* total number of grid points */
 28:   PetscReal   h;                     /* mesh width h = 1/(M-1) */
 29:   PetscReal   norm_2,norm_max;
 30:   PetscTruth  nox;                   /* indicates problem is to be run without graphics */
 31: } AppCtx;


 41: #define linear_no_matrix       0
 42: #define linear_no_time         1
 43: #define linear                 2
 44: #define nonlinear_no_jacobian  3
 45: #define nonlinear              4

 49: int main(int argc,char **argv)
 50: {
 52:   PetscInt       time_steps = 100,steps,m;
 53:   PetscMPIInt    size;
 54:   PetscInt       problem = linear_no_matrix;
 55:   PetscTruth     flg;
 56:   AppCtx         appctx;
 57:   PetscReal      dt,ftime;
 58:   TS             ts;
 59:   Mat            A=0,Alhs=0;
 60:   MatStructure   A_structure;
 61:   TSProblemType  tsproblem = TS_LINEAR;
 62:   PetscDraw      draw;
 63:   PetscViewer    viewer;
 64:   char           tsinfo[120];
 65:   TSType         type;
 66:   PetscTruth     isBEULER;
 67: 
 68:   PetscInitialize(&argc,&argv,(char*)0,help);
 69:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

 71:   appctx.M = 60;
 72:   PetscOptionsGetInt(PETSC_NULL,"-M",&appctx.M,PETSC_NULL);
 73:   PetscOptionsGetInt(PETSC_NULL,"-time",&time_steps,PETSC_NULL);
 74: 
 75:   PetscOptionsHasName(PETSC_NULL,"-nox",&appctx.nox);
 76:   appctx.norm_2 = 0.0; appctx.norm_max = 0.0;

 78:   /* Set up the ghost point communication pattern */
 79:   DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,appctx.M,1,1,PETSC_NULL,&appctx.da);
 80:   DACreateGlobalVector(appctx.da,&appctx.global);
 81:   VecGetLocalSize(appctx.global,&m);
 82:   DACreateLocalVector(appctx.da,&appctx.local);

 84:   /* Set up display to show wave graph */

 86:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,380,400,160,&appctx.viewer1);
 87:   PetscViewerDrawGetDraw(appctx.viewer1,0,&draw);
 88:   PetscDrawSetDoubleBuffer(draw);
 89:   PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,0,400,160,&appctx.viewer2);
 90:   PetscViewerDrawGetDraw(appctx.viewer2,0,&draw);
 91:   PetscDrawSetDoubleBuffer(draw);


 94:   /* make work array for evaluating right hand side function */
 95:   VecDuplicate(appctx.local,&appctx.localwork);

 97:   /* make work array for storing exact solution */
 98:   VecDuplicate(appctx.global,&appctx.solution);

100:   appctx.h = 1.0/(appctx.M-1.0);

102:   /* set initial conditions */
103:   Initial(appctx.global,&appctx);
104: 
105:   /*
106:      This example is written to allow one to easily test parts 
107:     of TS, we do not expect users to generally need to use more
108:     then a single TSProblemType
109:   */
110:   PetscOptionsHasName(PETSC_NULL,"-linear_no_matrix",&flg);
111:   if (flg) {
112:     tsproblem = TS_LINEAR;
113:     problem   = linear_no_matrix;
114:   }
115:   PetscOptionsHasName(PETSC_NULL,"-linear_constant_matrix",&flg);
116:   if (flg) {
117:     tsproblem = TS_LINEAR;
118:     problem   = linear_no_time;
119:   }
120:   PetscOptionsHasName(PETSC_NULL,"-linear_variable_matrix",&flg);
121:   if (flg) {
122:     tsproblem = TS_LINEAR;
123:     problem   = linear;
124:   }
125:   PetscOptionsHasName(PETSC_NULL,"-nonlinear_no_jacobian",&flg);
126:   if (flg) {
127:     tsproblem = TS_NONLINEAR;
128:     problem   = nonlinear_no_jacobian;
129:   }
130:   PetscOptionsHasName(PETSC_NULL,"-nonlinear_jacobian",&flg);
131:   if (flg) {
132:     tsproblem = TS_NONLINEAR;
133:     problem   = nonlinear;
134:   }
135: 
136:   /* make timestep context */
137:   TSCreate(PETSC_COMM_WORLD,&ts);
138:   TSSetProblemType(ts,tsproblem);
139:   TSSetMonitor(ts,Monitor,&appctx,PETSC_NULL);

141:   dt = appctx.h*appctx.h/2.01;

143:   if (problem == linear_no_matrix) {
144:     /*
145:          The user provides the RHS as a Shell matrix.
146:     */
147:     MatCreateShell(PETSC_COMM_WORLD,m,appctx.M,appctx.M,appctx.M,&appctx,&A);
148:     MatShellSetOperation(A,MATOP_MULT,(void(*)(void))RHSMatrixFree);
149:     TSSetRHSMatrix(ts,A,A,PETSC_NULL,&appctx);

151:   } else if (problem == linear_no_time) {
152:     /*
153:          The user provides the RHS as a constant matrix
154:     */
155:     MatCreate(PETSC_COMM_WORLD,&A);
156:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,appctx.M,appctx.M);
157:     MatSetFromOptions(A);
158:     RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
159:     TSSetRHSMatrix(ts,A,A,PETSC_NULL,&appctx);
160:   } else if (problem == linear) {
161:     /*
162:          The user provides the RHS as a time dependent matrix
163:     */
164:     MatCreate(PETSC_COMM_WORLD,&A);
165:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,appctx.M,appctx.M);
166:     MatSetFromOptions(A);
167:     RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
168:     TSSetRHSMatrix(ts,A,A,RHSMatrixHeat,&appctx);
169:   } else if (problem == nonlinear_no_jacobian) {
170:     /*
171:          The user provides the RHS and a Shell Jacobian
172:     */
173:     TSSetRHSFunction(ts,RHSFunctionHeat,&appctx);
174:     MatCreateShell(PETSC_COMM_WORLD,m,appctx.M,appctx.M,appctx.M,&appctx,&A);
175:     MatShellSetOperation(A,MATOP_MULT,(void(*)(void))RHSMatrixFree);
176:     TSSetRHSJacobian(ts,A,A,PETSC_NULL,&appctx);
177:   } else if (problem == nonlinear) {
178:     /*
179:          The user provides the RHS and Jacobian
180:     */
181:     TSSetRHSFunction(ts,RHSFunctionHeat,&appctx);
182:     MatCreate(PETSC_COMM_WORLD,&A);
183:     MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,appctx.M,appctx.M);
184:     MatSetFromOptions(A);
185:     RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
186:     TSSetRHSJacobian(ts,A,A,RHSJacobianHeat,&appctx);
187:   }

189:   TSSetFromOptions(ts);

191:   /* The user also provides the LHS matrix when ts_type = beuler */
192:   TSGetType(ts, &type);
193:   PetscStrcmp(type,"beuler",&isBEULER);
194:   if (isBEULER){
195:     if (problem == linear_no_time) {
196:       MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Alhs);
197:       MatZeroEntries(Alhs);
198:       MatShift(Alhs,1.0);
199:       TSSetLHSMatrix(ts,Alhs,Alhs,PETSC_NULL,PETSC_NULL);
200:     } else if (problem == linear) {
201:       MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Alhs);
202:       MatZeroEntries(Alhs);
203:       MatShift(Alhs,1.0);
204:       LHSMatrixHeat(ts,0.0,&Alhs,&Alhs,&A_structure,&appctx);
205:       TSSetLHSMatrix(ts,Alhs,Alhs,LHSMatrixHeat,&appctx);
206:     }
207:   }
208: 
209:   TSSetInitialTimeStep(ts,0.0,dt);
210:   TSSetDuration(ts,time_steps,100.);
211:   TSSetSolution(ts,appctx.global);


214:   TSSetUp(ts);
215:   TSStep(ts,&steps,&ftime);
216:   PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
217:   TSView(ts,viewer);

219:   PetscOptionsHasName(PETSC_NULL,"-test",&flg);
220:   if (flg) {
221:     PetscTruth iseuler;
222:     PetscTypeCompare((PetscObject)ts,"euler",&iseuler);
223:     if (iseuler) {
224:       if (!PETSC_NEAR(appctx.norm_2/steps,0.00257244,1.e-4)) {
225:         fprintf(stdout,"Error in Euler method: 2-norm %G expecting: 0.00257244\n",appctx.norm_2/steps);
226:       }
227:     } else {
228:       if (!PETSC_NEAR(appctx.norm_2/steps,0.00506174,1.e-4)) {
229:         fprintf(stdout,"Error in %s method: 2-norm %G expecting: 0.00506174\n",tsinfo,appctx.norm_2/steps);
230:       }
231:     }
232:   } else {
233:     PetscPrintf(PETSC_COMM_WORLD,"%D Procs Avg. error 2 norm %G max norm %G %s\n",
234:                 size,appctx.norm_2/steps,appctx.norm_max/steps,tsinfo);
235:   }

237:   PetscViewerDestroy(viewer);
238:   TSDestroy(ts);
239:   PetscViewerDestroy(appctx.viewer1);
240:   PetscViewerDestroy(appctx.viewer2);
241:   VecDestroy(appctx.localwork);
242:   VecDestroy(appctx.solution);
243:   VecDestroy(appctx.local);
244:   VecDestroy(appctx.global);
245:   DADestroy(appctx.da);
246:   if (A) {ierr= MatDestroy(A);}
247:   if (Alhs) {ierr= MatDestroy(Alhs);}

249:   PetscFinalize();
250:   return 0;
251: }

253: /* -------------------------------------------------------------------*/
256: PetscErrorCode Initial(Vec global,void *ctx)
257: {
258:   AppCtx         *appctx = (AppCtx*) ctx;
259:   PetscScalar    *localptr,h = appctx->h;
260:   PetscInt       i,mybase,myend;

263:   /* determine starting point of each processor */
264:   VecGetOwnershipRange(global,&mybase,&myend);

266:   /* Initialize the array */
267:   VecGetArray(global,&localptr);
268:   for (i=mybase; i<myend; i++) {
269:     localptr[i-mybase] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
270:   }
271:   VecRestoreArray(global,&localptr);
272:   return 0;
273: }

277: /*
278:        Exact solution 
279: */
280: PetscErrorCode Solution(PetscReal t,Vec solution,void *ctx)
281: {
282:   AppCtx *       appctx = (AppCtx*) ctx;
283:   PetscScalar    *localptr,h = appctx->h,ex1,ex2,sc1,sc2;
284:   PetscInt       i,mybase,myend;

287:   /* determine starting point of each processor */
288:   VecGetOwnershipRange(solution,&mybase,&myend);

290:   ex1 = exp(-36.*PETSC_PI*PETSC_PI*t);
291:   ex2 = exp(-4.*PETSC_PI*PETSC_PI*t);
292:   sc1 = PETSC_PI*6.*h;                 sc2 = PETSC_PI*2.*h;
293:   VecGetArray(solution,&localptr);
294:   for (i=mybase; i<myend; i++) {
295:     localptr[i-mybase] = PetscSinScalar(sc1*(PetscReal)i)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i)*ex2;
296:   }
297:   VecRestoreArray(solution,&localptr);
298:   return 0;
299: }

303: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal ltime,Vec global,void *ctx)
304: {
305:   AppCtx         *appctx = (AppCtx*) ctx;
307:   PetscReal      norm_2,norm_max;
308:   MPI_Comm       comm;

310:   PetscObjectGetComm((PetscObject)ts,&comm);

312:   VecView(global,appctx->viewer2);

314:   Solution(ltime,appctx->solution,ctx);
315:   VecAXPY(appctx->solution,-1.0,global);
316:   VecNorm(appctx->solution,NORM_2,&norm_2);
317:   norm_2 = sqrt(appctx->h)*norm_2;
318:   VecNorm(appctx->solution,NORM_MAX,&norm_max);

320:   if (!appctx->nox) {
321:     PetscPrintf(comm,"timestep %D time %G norm of error %G %G\n",step,ltime,norm_2,norm_max);
322:   }

324:   appctx->norm_2   += norm_2;
325:   appctx->norm_max += norm_max;

327:   VecView(appctx->solution,appctx->viewer1);

329:   return 0;
330: }

332: /* -----------------------------------------------------------------------*/
335: PetscErrorCode RHSMatrixFree(Mat mat,Vec x,Vec y)
336: {
337:   PetscErrorCode  ierr;
338:   void            *ctx;

340:   MatShellGetContext(mat,(void **)&ctx);
341:   RHSFunctionHeat(0,0.0,x,y,ctx);
342:   return 0;
343: }

347: PetscErrorCode RHSFunctionHeat(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
348: {
349:   AppCtx         *appctx = (AppCtx*) ctx;
350:   DA             da = appctx->da;
351:   Vec            local = appctx->local,localwork = appctx->localwork;
353:   PetscInt       i,localsize;
354:   PetscScalar    *copyptr,*localptr,sc;

356:   /*Extract local array */
357:   DAGlobalToLocalBegin(da,globalin,INSERT_VALUES,local);
358:   DAGlobalToLocalEnd(da,globalin,INSERT_VALUES,local);
359:   VecGetArray(local,&localptr);

361:   /* Extract work vector */
362:   VecGetArray(localwork,&copyptr);

364:   /* Update Locally - Make array of new values */
365:   /* Note: For the first and last entry I copy the value */
366:   /* if this is an interior node it is irrelevant */
367:   sc = 1.0/(appctx->h*appctx->h);
368:   VecGetLocalSize(local,&localsize);
369:   copyptr[0] = localptr[0];
370:   for (i=1; i<localsize-1; i++) {
371:     copyptr[i] = sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);
372:   }
373:   copyptr[localsize-1] = localptr[localsize-1];
374:   VecRestoreArray(local,&localptr);
375:   VecRestoreArray(localwork,&copyptr);

377:   /* Local to Global */
378:   DALocalToGlobal(da,localwork,INSERT_VALUES,globalout);
379:   return 0;
380: }

382: /* ---------------------------------------------------------------------*/
385: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
386: {
387:   Mat            A = *AA;
388:   AppCtx         *appctx = (AppCtx*) ctx;
390:   PetscInt       i,mstart,mend,idx[3];
391:   PetscMPIInt    size,rank;
392:   PetscScalar    v[3],stwo = -2./(appctx->h*appctx->h),sone = -.5*stwo;

394:   *str = SAME_NONZERO_PATTERN;
395: 
396:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
397:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

399:   MatGetOwnershipRange(A,&mstart,&mend);
400:   if (mstart == 0) {
401:     v[0] = 1.0;
402:     MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
403:     mstart++;
404:   }
405:   if (mend == appctx->M) {
406:     mend--;
407:     v[0] = 1.0;
408:     MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
409:   }

411:   /*
412:      Construct matrice one row at a time
413:   */
414:   v[0] = sone; v[1] = stwo; v[2] = sone;
415:   for (i=mstart; i<mend; i++) {
416:     idx[0] = i-1; idx[1] = i; idx[2] = i+1;
417:     MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
418:   }

420:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
421:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
422:   return 0;
423: }

427: PetscErrorCode RHSJacobianHeat(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
428: {
429:   return RHSMatrixHeat(ts,t,AA,BB,str,ctx);
430: }

434: PetscErrorCode LHSMatrixHeat(TS ts,PetscReal t,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
435: {
436:   *str = SAME_NONZERO_PATTERN;
437:   /* printf("  LHSMatrixHeat t: %g\n",t); */
438:   return 0;
439: }