Actual source code: ex1.c
petsc-3.3-p7 2013-05-11
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 <petscdmda.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: DM 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: PetscBool nox; /* indicates problem is to be run without graphics */
30: } AppCtx;
32: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void *);
33: extern PetscErrorCode RHSFunctionHeat(TS,PetscReal,Vec,Vec,void*);
34: extern PetscErrorCode RHSMatrixFree(Mat,Vec,Vec);
35: extern PetscErrorCode Initial(Vec,void*);
36: extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Mat *,Mat *,MatStructure *,void *);
37: extern PetscErrorCode LHSMatrixHeat(TS,PetscReal,Mat *,Mat *,MatStructure *,void *);
38: extern PetscErrorCode RHSJacobianHeat(TS,PetscReal,Vec,Mat*,Mat*,MatStructure *,void*);
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: PetscBool 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: char tsinfo[120];
63:
64: PetscInitialize(&argc,&argv,(char*)0,help);
65: MPI_Comm_size(PETSC_COMM_WORLD,&size);
67: appctx.M = 60;
68: PetscOptionsGetInt(PETSC_NULL,"-M",&appctx.M,PETSC_NULL);
69: PetscOptionsGetInt(PETSC_NULL,"-time",&maxsteps,PETSC_NULL);
70:
71: PetscOptionsHasName(PETSC_NULL,"-nox",&appctx.nox);
72: appctx.norm_2 = 0.0; appctx.norm_max = 0.0;
74: /* Set up the ghost point communication pattern */
75: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_NONE,appctx.M,1,1,PETSC_NULL,&appctx.da);
76: DMCreateGlobalVector(appctx.da,&appctx.global);
77: VecGetLocalSize(appctx.global,&m);
78: DMCreateLocalVector(appctx.da,&appctx.local);
80: /* Set up display to show wave graph */
81: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,380,400,160,&appctx.viewer1);
82: PetscViewerDrawGetDraw(appctx.viewer1,0,&draw);
83: PetscDrawSetDoubleBuffer(draw);
84: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,"",80,0,400,160,&appctx.viewer2);
85: PetscViewerDrawGetDraw(appctx.viewer2,0,&draw);
86: PetscDrawSetDoubleBuffer(draw);
88: /* make work array for evaluating right hand side function */
89: VecDuplicate(appctx.local,&appctx.localwork);
91: /* make work array for storing exact solution */
92: VecDuplicate(appctx.global,&appctx.solution);
94: appctx.h = 1.0/(appctx.M-1.0);
96: /* set initial conditions */
97: Initial(appctx.global,&appctx);
98:
99: /*
100: This example is written to allow one to easily test parts
101: of TS, we do not expect users to generally need to use more
102: then a single TSProblemType
103: */
104: tsproblem = TS_NONLINEAR;
105: problem = nonlinear_no_jacobian;
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_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: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This case needs rewritten");
144: /*
145: The user provides the RHS as a Shell matrix.
146: */
147: /*
148: MatCreateShell(PETSC_COMM_WORLD,m,m,PETSC_DECIDE,PETSC_DECIDE,&appctx,&A);
149: MatShellSetOperation(A,MATOP_MULT,(void(*)(void))RHSMatrixFree);
150: TSSetMatrices(ts,A,PETSC_NULL,PETSC_NULL,PETSC_NULL,DIFFERENT_NONZERO_PATTERN,&appctx);
151: */
152: } else if (problem == linear_constant_matrix) {
153: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This case needs rewritten");
154: /*
155: The user provides the RHS as a constant matrix
156: */
157: /*
158: MatCreate(PETSC_COMM_WORLD,&A);
159: MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE);
160: MatSetFromOptions(A);
161: RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
163: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Alhs);
164: MatZeroEntries(Alhs);
165: MatShift(Alhs,1.0);
166: TSSetMatrices(ts,A,PETSC_NULL,Alhs,PETSC_NULL,DIFFERENT_NONZERO_PATTERN,&appctx);
167: */
168: } else if (problem == linear_variable_matrix) {
169: SETERRQ(PETSC_COMM_SELF,PETSC_ERR_SUP,"This case needs rewritten");
170: /*
171: The user provides the RHS as a time dependent matrix
172: */
173: /*
174: MatCreate(PETSC_COMM_WORLD,&A);
175: MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE);
176: MatSetFromOptions(A);
177: RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
179: MatDuplicate(A,MAT_DO_NOT_COPY_VALUES,&Alhs);
180: MatZeroEntries(Alhs);
181: MatShift(Alhs,1.0);
182: LHSMatrixHeat(ts,0.0,&Alhs,&Alhs,&A_structure,&appctx);
183: TSSetMatrices(ts,A,RHSMatrixHeat,Alhs,LHSMatrixHeat,DIFFERENT_NONZERO_PATTERN,&appctx);
184: */
185: } else if (problem == nonlinear_no_jacobian) {
186: /*
187: The user provides the RHS and a Shell Jacobian
188: */
189: TSSetRHSFunction(ts,PETSC_NULL,RHSFunctionHeat,&appctx);
190: MatCreateShell(PETSC_COMM_WORLD,m,m,PETSC_DECIDE,PETSC_DECIDE,&appctx,&A);
191: MatShellSetOperation(A,MATOP_MULT,(void(*)(void))RHSMatrixFree);
192: TSSetRHSJacobian(ts,A,A,PETSC_NULL,&appctx);
193: } else if (problem == nonlinear_jacobian) {
194: /*
195: The user provides the RHS and Jacobian
196: */
197: TSSetRHSFunction(ts,PETSC_NULL,RHSFunctionHeat,&appctx);
198: MatCreate(PETSC_COMM_WORLD,&A);
199: MatSetSizes(A,m,m,PETSC_DECIDE,PETSC_DECIDE);
200: MatSetFromOptions(A);
201: MatSetUp(A);
202: RHSMatrixHeat(ts,0.0,&A,&A,&A_structure,&appctx);
203: TSSetRHSJacobian(ts,A,A,RHSJacobianHeat,&appctx);
204: }
206: TSSetInitialTimeStep(ts,0.0,dt);
207: TSSetDuration(ts,maxsteps,maxtime);
208: TSSetSolution(ts,appctx.global);
210: TSSetFromOptions(ts);
212: TSSolve(ts,appctx.global,&ftime);
214: PetscOptionsHasName(PETSC_NULL,"-testinfo",&flg);
215: if (flg) {
216: PetscBool iseuler;
217: TSGetTimeStepNumber(ts,&steps);
219: PetscObjectTypeCompare((PetscObject)ts,"euler",&iseuler);
220: if (iseuler) {
221: if (!PETSC_NEAR(appctx.norm_2/steps,0.00257244,1.e-4)) {
222: fprintf(stdout,"Error in Euler method: 2-norm %G expecting: 0.00257244\n",appctx.norm_2/steps);
223: }
224: } else {
225: PetscPrintf(PETSC_COMM_WORLD,"%D Procs; Avg. error 2-norm %G; max-norm %G; %s\n",
226: size,appctx.norm_2/steps,appctx.norm_max/steps,tsinfo);
227: }
228: }
230: TSDestroy(&ts);
231: PetscViewerDestroy(&appctx.viewer1);
232: PetscViewerDestroy(&appctx.viewer2);
233: VecDestroy(&appctx.localwork);
234: VecDestroy(&appctx.solution);
235: VecDestroy(&appctx.local);
236: VecDestroy(&appctx.global);
237: DMDestroy(&appctx.da);
238: if (A) {ierr= MatDestroy(&A);}
239: if (Alhs) {ierr= MatDestroy(&Alhs);}
241: PetscFinalize();
242: return 0;
243: }
245: /* -------------------------------------------------------------------*/
246: /*
247: Set initial condition: u(t=0) = sin(6*pi*x) + 3*sin(2*pi*x)
248: */
251: PetscErrorCode Initial(Vec global,void *ctx)
252: {
253: AppCtx *appctx = (AppCtx*) ctx;
254: PetscScalar *localptr,h = appctx->h;
255: PetscInt i,mybase,myend;
259: /* determine starting point of each processor */
260: VecGetOwnershipRange(global,&mybase,&myend);
262: /* Initialize the array */
263: VecGetArray(global,&localptr);
264: for (i=mybase; i<myend; i++) {
265: localptr[i-mybase] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
266: }
267: VecRestoreArray(global,&localptr);
268: return(0);
269: }
273: /*
274: Exact solution:
275: u = sin(6*pi*x)*exp(-36*pi*pi*t) + 3*sin(2*pi*x)*exp(-4*pi*pi*t)
276: */
277: PetscErrorCode Solution(PetscReal t,Vec solution,void *ctx)
278: {
279: AppCtx * appctx = (AppCtx*) ctx;
280: PetscScalar *localptr,h = appctx->h,ex1,ex2,sc1,sc2;
281: PetscInt i,mybase,myend;
285: /* determine starting point of each processor */
286: VecGetOwnershipRange(solution,&mybase,&myend);
288: ex1 = exp(-36.*PETSC_PI*PETSC_PI*t);
289: ex2 = exp(-4.*PETSC_PI*PETSC_PI*t);
290: sc1 = PETSC_PI*6.*h;
291: sc2 = PETSC_PI*2.*h;
292: VecGetArray(solution,&localptr);
293: for (i=mybase; i<myend; i++) {
294: localptr[i-mybase] = PetscSinScalar(sc1*(PetscReal)i)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i)*ex2;
295: }
296: VecRestoreArray(solution,&localptr);
297: return(0);
298: }
300: /*
301: step - iteration number
302: ltime - current time
303: global - current iterate
304: */
307: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal ltime,Vec global,void *ctx)
308: {
309: AppCtx *appctx = (AppCtx*) ctx;
311: PetscReal norm_2,norm_max;
312: MPI_Comm comm;
315: PetscObjectGetComm((PetscObject)ts,&comm);
316: if (!appctx->nox) {
317: VecView(global,appctx->viewer2); /* show wave graph */
318: }
319: Solution(ltime,appctx->solution,ctx); /* get true solution at current time */
320: VecAXPY(appctx->solution,-1.0,global);
321: VecNorm(appctx->solution,NORM_2,&norm_2);
322: norm_2 = PetscSqrtReal(appctx->h)*norm_2;
323: VecNorm(appctx->solution,NORM_MAX,&norm_max);
324: PetscPrintf(comm,"timestep %D time %G norm of error %.5f %.5f\n",step,ltime,norm_2,norm_max);
326: appctx->norm_2 += norm_2;
327: appctx->norm_max += norm_max;
328: if (!appctx->nox) {
329: VecView(appctx->solution,appctx->viewer1);
330: }
331: return(0);
332: }
334: /* -----------------------------------------------------------------------*/
337: PetscErrorCode RHSMatrixFree(Mat mat,Vec x,Vec y)
338: {
339: PetscErrorCode ierr;
340: void *ctx;
343: MatShellGetContext(mat,(void **)&ctx);
344: RHSFunctionHeat(0,0.0,x,y,ctx);
345: return(0);
346: }
350: PetscErrorCode RHSFunctionHeat(TS ts,PetscReal t,Vec globalin,Vec globalout,void *ctx)
351: {
352: AppCtx *appctx = (AppCtx*) ctx;
353: DM da = appctx->da;
354: Vec local = appctx->local,localwork = appctx->localwork;
356: PetscInt i,localsize;
357: PetscScalar *copyptr,*localptr,sc;
360: /*Extract local array */
361: DMGlobalToLocalBegin(da,globalin,INSERT_VALUES,local);
362: DMGlobalToLocalEnd(da,globalin,INSERT_VALUES,local);
363: VecGetArray(local,&localptr);
365: /* Extract work vector */
366: VecGetArray(localwork,©ptr);
368: /* Update Locally - Make array of new values */
369: /* Note: For the first and last entry I copy the value */
370: /* if this is an interior node it is irrelevant */
371: sc = 1.0/(appctx->h*appctx->h);
372: VecGetLocalSize(local,&localsize);
373: copyptr[0] = localptr[0];
374: for (i=1; i<localsize-1; i++) {
375: copyptr[i] = sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);
376: }
377: copyptr[localsize-1] = localptr[localsize-1];
378: VecRestoreArray(local,&localptr);
379: VecRestoreArray(localwork,©ptr);
381: /* Local to Global */
382: DMLocalToGlobalBegin(da,localwork,INSERT_VALUES,globalout);
383: DMLocalToGlobalEnd(da,localwork,INSERT_VALUES,globalout);
384: return(0);
385: }
387: /* ---------------------------------------------------------------------*/
390: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
391: {
392: Mat A = *AA;
393: AppCtx *appctx = (AppCtx*) ctx;
395: PetscInt i,mstart,mend,idx[3];
396: PetscMPIInt size,rank;
397: PetscScalar v[3],stwo = -2./(appctx->h*appctx->h),sone = -.5*stwo;
400: *str = SAME_NONZERO_PATTERN;
401:
402: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
403: MPI_Comm_size(PETSC_COMM_WORLD,&size);
405: MatGetOwnershipRange(A,&mstart,&mend);
406: if (mstart == 0) {
407: v[0] = 1.0;
408: MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
409: mstart++;
410: }
411: if (mend == appctx->M) {
412: mend--;
413: v[0] = 1.0;
414: MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
415: }
417: /*
418: Construct matrice one row at a time
419: */
420: v[0] = sone; v[1] = stwo; v[2] = sone;
421: for (i=mstart; i<mend; i++) {
422: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
423: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
424: }
426: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
427: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
428: return(0);
429: }
433: PetscErrorCode RHSJacobianHeat(TS ts,PetscReal t,Vec x,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
434: {
436: RHSMatrixHeat(ts,t,AA,BB,str,ctx);
437: return(0);
438: }
440: /* A = indentity matrix */
443: PetscErrorCode LHSMatrixHeat(TS ts,PetscReal t,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
444: {
446: Mat A = *AA;
449: *str = SAME_NONZERO_PATTERN;
450: MatZeroEntries(A);
451: MatShift(A,1.0);
452: return(0);
453: }