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 -pc_type lu
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_constant_matrix 1
43: #define linear_variable_matrix 2
44: #define nonlinear_no_jacobian 3
45: #define nonlinear_jacobian 4
49: int main(int argc,char **argv)
50: {
52: PetscInt maxsteps = 100,steps,m;
53: PetscMPIInt size;
54: PetscInt problem = linear_no_matrix;
55: PetscTruth flg;
56: AppCtx appctx;
57: PetscReal dt,ftime,maxtime=100.;
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:
66: PetscInitialize(&argc,&argv,(char*)0,help);
67: MPI_Comm_size(PETSC_COMM_WORLD,&size);
69: appctx.M = 60;
70: PetscOptionsGetInt(PETSC_NULL,"-M",&appctx.M,PETSC_NULL);
71: PetscOptionsGetInt(PETSC_NULL,"-time",&maxsteps,PETSC_NULL);
72:
73: PetscOptionsHasName(PETSC_NULL,"-nox",&appctx.nox);
74: appctx.norm_2 = 0.0; appctx.norm_max = 0.0;
76: /* Set up the ghost point communication pattern */
77: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,appctx.M,1,1,PETSC_NULL,&appctx.da);
78: DACreateGlobalVector(appctx.da,&appctx.global);
79: VecGetLocalSize(appctx.global,&m);
80: DACreateLocalVector(appctx.da,&appctx.local);
82: /* Set up display to show wave graph */
83: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,380,400,160,&appctx.viewer1);
84: PetscViewerDrawGetDraw(appctx.viewer1,0,&draw);
85: PetscDrawSetDoubleBuffer(draw);
86: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,0,400,160,&appctx.viewer2);
87: PetscViewerDrawGetDraw(appctx.viewer2,0,&draw);
88: PetscDrawSetDoubleBuffer(draw);
90: /* make work array for evaluating right hand side function */
91: VecDuplicate(appctx.local,&appctx.localwork);
93: /* make work array for storing exact solution */
94: VecDuplicate(appctx.global,&appctx.solution);
96: appctx.h = 1.0/(appctx.M-1.0);
98: /* set initial conditions */
99: Initial(appctx.global,&appctx);
100:
101: /*
102: This example is written to allow one to easily test parts
103: of TS, we do not expect users to generally need to use more
104: then a single TSProblemType
105: */
106: PetscOptionsHasName(PETSC_NULL,"-linear_no_matrix",&flg);
107: if (flg) {
108: tsproblem = TS_LINEAR;
109: problem = linear_no_matrix;
110: }
111: PetscOptionsHasName(PETSC_NULL,"-linear_constant_matrix",&flg);
112: if (flg) {
113: tsproblem = TS_LINEAR;
114: problem = linear_constant_matrix;
115: }
116: PetscOptionsHasName(PETSC_NULL,"-linear_variable_matrix",&flg);
117: if (flg) {
118: tsproblem = TS_LINEAR;
119: problem = linear_variable_matrix;
120: }
121: PetscOptionsHasName(PETSC_NULL,"-nonlinear_no_jacobian",&flg);
122: if (flg) {
123: tsproblem = TS_NONLINEAR;
124: problem = nonlinear_no_jacobian;
125: }
126: PetscOptionsHasName(PETSC_NULL,"-nonlinear_jacobian",&flg);
127: if (flg) {
128: tsproblem = TS_NONLINEAR;
129: problem = nonlinear_jacobian;
130: }
131:
132: /* make timestep context */
133: TSCreate(PETSC_COMM_WORLD,&ts);
134: TSSetProblemType(ts,tsproblem);
135: PetscOptionsHasName(PETSC_NULL,"-monitor",&flg);
136: if (flg) {
137: TSMonitorSet(ts,Monitor,&appctx,PETSC_NULL);
138: }
140: dt = appctx.h*appctx.h/2.01;
142: if (problem == linear_no_matrix) {
143: /*
144: The user provides the RHS as a Shell matrix.
145: */
146: MatCreateShell(PETSC_COMM_WORLD,m,m,PETSC_DECIDE,PETSC_DECIDE,&appctx,&A);
147: MatShellSetOperation(A,MATOP_MULT,(void(*)(void))RHSMatrixFree);
148: TSSetMatrices(ts,A,PETSC_NULL,PETSC_NULL,PETSC_NULL,DIFFERENT_NONZERO_PATTERN,&appctx);
149: } else if (problem == linear_constant_matrix) {
150: /*
151: The user provides the RHS as a constant matrix
152: */
153: MatCreate(PETSC_COMM_WORLD,&A);
154: MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE);
155: MatSetFromOptions(A);
156: RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx); /* A is assembled here */
158: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Alhs);
159: MatZeroEntries(Alhs);
160: MatShift(Alhs,1.0);
161: TSSetMatrices(ts,A,PETSC_NULL,Alhs,PETSC_NULL,DIFFERENT_NONZERO_PATTERN,&appctx);
162: } else if (problem == linear_variable_matrix) {
163: /*
164: The user provides the RHS as a time dependent matrix
165: */
166: MatCreate(PETSC_COMM_WORLD,&A);
167: MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE);
168: MatSetFromOptions(A);
169: RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
171: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Alhs);
172: MatZeroEntries(Alhs);
173: MatShift(Alhs,1.0);
174: LHSMatrixHeat(ts,0.0,&Alhs,&Alhs,&A_structure,&appctx);
175: TSSetMatrices(ts,A,RHSMatrixHeat,Alhs,LHSMatrixHeat,DIFFERENT_NONZERO_PATTERN,&appctx);
176: } else if (problem == nonlinear_no_jacobian) {
177: /*
178: The user provides the RHS and a Shell Jacobian
179: */
180: TSSetRHSFunction(ts,RHSFunctionHeat,&appctx);
181: MatCreateShell(PETSC_COMM_WORLD,m,m,PETSC_DECIDE,PETSC_DECIDE,&appctx,&A);
182: MatShellSetOperation(A,MATOP_MULT,(void(*)(void))RHSMatrixFree);
183: TSSetRHSJacobian(ts,A,A,PETSC_NULL,&appctx);
184: } else if (problem == nonlinear_jacobian) {
185: /*
186: The user provides the RHS and Jacobian
187: */
188: TSSetRHSFunction(ts,RHSFunctionHeat,&appctx);
189: MatCreate(PETSC_COMM_WORLD,&A);
190: MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE);
191: MatSetFromOptions(A);
192: RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
193: TSSetRHSJacobian(ts,A,A,RHSJacobianHeat,&appctx);
194: }
195: TSSetFromOptions(ts);
197: TSSetInitialTimeStep(ts,0.0,dt);
198: TSSetDuration(ts,maxsteps,maxtime);
199: TSSetSolution(ts,appctx.global);
201: TSSetUp(ts);
202: TSStep(ts,&steps,&ftime);
203: PetscViewerStringOpen(PETSC_COMM_WORLD,tsinfo,120,&viewer);
204: if (size == 1){ /* TSView() crashes for non euler methods with np>1 ? */
205: TSView(ts,viewer);
206: }
208: PetscOptionsHasName(PETSC_NULL,"-testinfo",&flg);
209: if (flg) {
210: PetscTruth iseuler;
211: PetscTypeCompare((PetscObject)ts,"euler",&iseuler);
212: if (iseuler) {
213: if (!PETSC_NEAR(appctx.norm_2/steps,0.00257244,1.e-4)) {
214: fprintf(stdout,"Error in Euler method: 2-norm %G expecting: 0.00257244\n",appctx.norm_2/steps);
215: }
216: } else {
217: PetscPrintf(PETSC_COMM_WORLD,"%D Procs; Avg. error 2-norm %G; max-norm %G; %s\n",
218: size,appctx.norm_2/steps,appctx.norm_max/steps,tsinfo);
219: }
220: }
222: PetscViewerDestroy(viewer);
223: TSDestroy(ts);
224: PetscViewerDestroy(appctx.viewer1);
225: PetscViewerDestroy(appctx.viewer2);
226: VecDestroy(appctx.localwork);
227: VecDestroy(appctx.solution);
228: VecDestroy(appctx.local);
229: VecDestroy(appctx.global);
230: DADestroy(appctx.da);
231: if (A) {ierr= MatDestroy(A);}
232: if (Alhs) {ierr= MatDestroy(Alhs);}
234: PetscFinalize();
235: return 0;
236: }
238: /* -------------------------------------------------------------------*/
239: /*
240: Set initial condition: u(t=0) = sin(6*pi*x) + 3*sin(2*pi*x)
241: */
244: PetscErrorCode Initial(Vec global,void *ctx)
245: {
246: AppCtx *appctx = (AppCtx*) ctx;
247: PetscScalar *localptr,h = appctx->h;
248: PetscInt i,mybase,myend;
252: /* determine starting point of each processor */
253: VecGetOwnershipRange(global,&mybase,&myend);
255: /* Initialize the array */
256: VecGetArray(global,&localptr);
257: for (i=mybase; i<myend; i++) {
258: localptr[i-mybase] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
259: }
260: VecRestoreArray(global,&localptr);
261: return(0);
262: }
266: /*
267: Exact solution:
268: u = sin(6*pi*x)*exp(-36*pi*pi*t) + 3*sin(2*pi*x)*exp(-4*pi*pi*t)
269: */
270: PetscErrorCode Solution(PetscReal t,Vec solution,void *ctx)
271: {
272: AppCtx * appctx = (AppCtx*) ctx;
273: PetscScalar *localptr,h = appctx->h,ex1,ex2,sc1,sc2;
274: PetscInt i,mybase,myend;
278: /* determine starting point of each processor */
279: VecGetOwnershipRange(solution,&mybase,&myend);
281: ex1 = exp(-36.*PETSC_PI*PETSC_PI*t);
282: ex2 = exp(-4.*PETSC_PI*PETSC_PI*t);
283: sc1 = PETSC_PI*6.*h;
284: sc2 = PETSC_PI*2.*h;
285: VecGetArray(solution,&localptr);
286: for (i=mybase; i<myend; i++) {
287: localptr[i-mybase] = PetscSinScalar(sc1*(PetscReal)i)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i)*ex2;
288: }
289: VecRestoreArray(solution,&localptr);
290: return(0);
291: }
293: /*
294: step - iteration number
295: ltime - current time
296: global - current iterate
297: */
300: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal ltime,Vec global,void *ctx)
301: {
302: AppCtx *appctx = (AppCtx*) ctx;
304: PetscReal norm_2,norm_max;
305: MPI_Comm comm;
308: PetscObjectGetComm((PetscObject)ts,&comm);
309: if (!appctx->nox) {
310: VecView(global,appctx->viewer2); /* show wave graph */
311: }
312: Solution(ltime,appctx->solution,ctx); /* get true solution at current time */
313: VecAXPY(appctx->solution,-1.0,global);
314: VecNorm(appctx->solution,NORM_2,&norm_2);
315: norm_2 = sqrt(appctx->h)*norm_2;
316: VecNorm(appctx->solution,NORM_MAX,&norm_max);
317: PetscPrintf(comm,"timestep %D time %G norm of error %.5f %.5f\n",step,ltime,norm_2,norm_max);
319: appctx->norm_2 += norm_2;
320: appctx->norm_max += norm_max;
321: if (!appctx->nox) {
322: VecView(appctx->solution,appctx->viewer1);
323: }
324: return(0);
325: }
327: /* -----------------------------------------------------------------------*/
330: PetscErrorCode RHSMatrixFree(Mat mat,Vec x,Vec y)
331: {
332: PetscErrorCode ierr;
333: void *ctx;
336: MatShellGetContext(mat,(void **)&ctx);
337: RHSFunctionHeat(0,0.0,x,y,ctx);
338: return(0);
339: }
343: PetscErrorCode RHSFunctionHeat(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
344: {
345: AppCtx *appctx = (AppCtx*) ctx;
346: DA da = appctx->da;
347: Vec local = appctx->local,localwork = appctx->localwork;
349: PetscInt i,localsize;
350: PetscScalar *copyptr,*localptr,sc;
353: /*Extract local array */
354: DAGlobalToLocalBegin(da,globalin,INSERT_VALUES,local);
355: DAGlobalToLocalEnd(da,globalin,INSERT_VALUES,local);
356: VecGetArray(local,&localptr);
358: /* Extract work vector */
359: VecGetArray(localwork,©ptr);
361: /* Update Locally - Make array of new values */
362: /* Note: For the first and last entry I copy the value */
363: /* if this is an interior node it is irrelevant */
364: sc = 1.0/(appctx->h*appctx->h);
365: VecGetLocalSize(local,&localsize);
366: copyptr[0] = localptr[0];
367: for (i=1; i<localsize-1; i++) {
368: copyptr[i] = sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);
369: }
370: copyptr[localsize-1] = localptr[localsize-1];
371: VecRestoreArray(local,&localptr);
372: VecRestoreArray(localwork,©ptr);
374: /* Local to Global */
375: DALocalToGlobal(da,localwork,INSERT_VALUES,globalout);
376: return(0);
377: }
379: /* ---------------------------------------------------------------------*/
382: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
383: {
384: Mat A = *AA;
385: AppCtx *appctx = (AppCtx*) ctx;
387: PetscInt i,mstart,mend,idx[3];
388: PetscMPIInt size,rank;
389: PetscScalar v[3],stwo = -2./(appctx->h*appctx->h),sone = -.5*stwo;
392: *str = SAME_NONZERO_PATTERN;
393:
394: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
395: MPI_Comm_size(PETSC_COMM_WORLD,&size);
397: MatGetOwnershipRange(A,&mstart,&mend);
398: if (mstart == 0) {
399: v[0] = 1.0;
400: MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
401: mstart++;
402: }
403: if (mend == appctx->M) {
404: mend--;
405: v[0] = 1.0;
406: MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
407: }
409: /*
410: Construct matrice one row at a time
411: */
412: v[0] = sone; v[1] = stwo; v[2] = sone;
413: for (i=mstart; i<mend; i++) {
414: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
415: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
416: }
418: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
419: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
420: return(0);
421: }
425: PetscErrorCode RHSJacobianHeat(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
426: {
428: RHSMatrixHeat(ts,t,AA,BB,str,ctx);
429: return(0);
430: }
432: /* A = indentity matrix */
435: PetscErrorCode LHSMatrixHeat(TS ts,PetscReal t,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
436: {
438: Mat A = *AA;
441: *str = SAME_NONZERO_PATTERN;
442: MatZeroEntries(A);
443: MatShift(A,1.0);
444: return(0);
445: }