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