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,©ptr);
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,©ptr);
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: }