Actual source code: ex2.c
2: static char help[] ="Solves a time-dependent nonlinear PDE. Uses implicit\n\
3: timestepping. Runtime options include:\n\
4: -M <xg>, where <xg> = number of grid points\n\
5: -debug : Activate debugging printouts\n\
6: -nox : Deactivate x-window graphics\n\n";
8: /*
9: Concepts: TS^time-dependent nonlinear problems
10: Processors: n
11: */
13: /* ------------------------------------------------------------------------
15: This program solves the PDE
17: u * u_xx
18: u_t = ---------
19: 2*(t+1)^2
21: on the domain 0 <= x <= 1, with boundary conditions
22: u(t,0) = t + 1, u(t,1) = 2*t + 2,
23: and initial condition
24: u(0,x) = 1 + x*x.
26: The exact solution is:
27: u(t,x) = (1 + x*x) * (1 + t)
29: Note that since the solution is linear in time and quadratic in x,
30: the finite difference scheme actually computes the "exact" solution.
32: We use by default the backward Euler method.
34: ------------------------------------------------------------------------- */
36: /*
37: Include "petscts.h" to use the PETSc timestepping routines. Note that
38: this file automatically includes "petsc.h" and other lower-level
39: PETSc include files.
41: Include the "petscda.h" to allow us to use the distributed array data
42: structures to manage the parallel grid.
43: */
44: #include petscts.h
45: #include petscda.h
47: /*
48: User-defined application context - contains data needed by the
49: application-provided callback routines.
50: */
51: typedef struct {
52: MPI_Comm comm; /* communicator */
53: DA da; /* distributed array data structure */
54: Vec localwork; /* local ghosted work vector */
55: Vec u_local; /* local ghosted approximate solution vector */
56: Vec solution; /* global exact solution vector */
57: PetscInt m; /* total number of grid points */
58: PetscReal h; /* mesh width: h = 1/(m-1) */
59: PetscTruth debug; /* flag (1 indicates activation of debugging printouts) */
60: } AppCtx;
62: /*
63: User-defined routines, provided below.
64: */
71: /*
72: Utility routine for finite difference Jacobian approximation
73: */
78: int main(int argc,char **argv)
79: {
80: AppCtx appctx; /* user-defined application context */
81: TS ts; /* timestepping context */
82: Mat A; /* Jacobian matrix data structure */
83: Vec u; /* approximate solution vector */
84: PetscInt time_steps_max = 1000; /* default max timesteps */
86: PetscInt steps;
87: PetscReal ftime; /* final time */
88: PetscReal dt;
89: PetscReal time_total_max = 100.0; /* default max total time */
90: PetscTruth flg;
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Initialize program and set problem parameters
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
95:
96: PetscInitialize(&argc,&argv,(char*)0,help);
98: appctx.comm = PETSC_COMM_WORLD;
99: appctx.m = 60;
100: PetscOptionsGetInt(PETSC_NULL,"-M",&appctx.m,PETSC_NULL);
101: PetscOptionsHasName(PETSC_NULL,"-debug",&appctx.debug);
102: appctx.h = 1.0/(appctx.m-1.0);
104: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
105: Create vector data structures
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: /*
109: Create distributed array (DA) to manage parallel grid and vectors
110: and to set up the ghost point communication pattern. There are M
111: total grid values spread equally among all the processors.
112: */
113: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,appctx.m,1,1,PETSC_NULL,
114: &appctx.da);
116: /*
117: Extract global and local vectors from DA; we use these to store the
118: approximate solution. Then duplicate these for remaining vectors that
119: have the same types.
120: */
121: DACreateGlobalVector(appctx.da,&u);
122: DACreateLocalVector(appctx.da,&appctx.u_local);
124: /*
125: Create local work vector for use in evaluating right-hand-side function;
126: create global work vector for storing exact solution.
127: */
128: VecDuplicate(appctx.u_local,&appctx.localwork);
129: VecDuplicate(u,&appctx.solution);
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: Create timestepping solver context; set callback routine for
133: right-hand-side function evaluation.
134: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
136: TSCreate(PETSC_COMM_WORLD,&ts);
137: TSSetProblemType(ts,TS_NONLINEAR);
138: TSSetRHSFunction(ts,RHSFunction,&appctx);
140: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
141: Set optional user-defined monitoring routine
142: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
144: TSMonitorSet(ts,Monitor,&appctx,PETSC_NULL);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: For nonlinear problems, the user can provide a Jacobian evaluation
148: routine (or use a finite differencing approximation).
150: Create matrix data structure; set Jacobian evaluation routine.
151: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
153: MatCreate(PETSC_COMM_WORLD,&A);
154: MatSetSizes(A,PETSC_DECIDE,PETSC_DECIDE,appctx.m,appctx.m);
155: MatSetFromOptions(A);
156: PetscOptionsHasName(PETSC_NULL,"-fdjac",&flg);
157: if (flg) {
158: TSSetRHSJacobian(ts,A,A,RHSJacobianFD,&appctx);
159: } else {
160: TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx);
161: }
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Set solution vector and initial timestep
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
167: dt = appctx.h/2.0;
168: TSSetInitialTimeStep(ts,0.0,dt);
169: TSSetSolution(ts,u);
171: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
172: Customize timestepping solver:
173: - Set the solution method to be the Backward Euler method.
174: - Set timestepping duration info
175: Then set runtime options, which can override these defaults.
176: For example,
177: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
178: to override the defaults set by TSSetDuration().
179: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: TSSetType(ts,TS_BEULER);
182: TSSetDuration(ts,time_steps_max,time_total_max);
183: TSSetFromOptions(ts);
185: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
186: Solve the problem
187: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
189: /*
190: Evaluate initial conditions
191: */
192: InitialConditions(u,&appctx);
194: /*
195: Run the timestepping solver
196: */
197: TSStep(ts,&steps,&ftime);
199: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
200: Free work space. All PETSc objects should be destroyed when they
201: are no longer needed.
202: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: TSDestroy(ts);
205: VecDestroy(u);
206: MatDestroy(A);
207: DADestroy(appctx.da);
208: VecDestroy(appctx.localwork);
209: VecDestroy(appctx.solution);
210: VecDestroy(appctx.u_local);
212: /*
213: Always call PetscFinalize() before exiting a program. This routine
214: - finalizes the PETSc libraries as well as MPI
215: - provides summary and diagnostic information if certain runtime
216: options are chosen (e.g., -log_summary).
217: */
218: PetscFinalize();
219: return 0;
220: }
221: /* --------------------------------------------------------------------- */
224: /*
225: InitialConditions - Computes the solution at the initial time.
227: Input Parameters:
228: u - uninitialized solution vector (global)
229: appctx - user-defined application context
231: Output Parameter:
232: u - vector with solution at initial time (global)
233: */
234: PetscErrorCode InitialConditions(Vec u,AppCtx *appctx)
235: {
236: PetscScalar *u_localptr,h = appctx->h,x;
237: PetscInt i,mybase,myend;
240: /*
241: Determine starting point of each processor's range of
242: grid values.
243: */
244: VecGetOwnershipRange(u,&mybase,&myend);
246: /*
247: Get a pointer to vector data.
248: - For default PETSc vectors, VecGetArray() returns a pointer to
249: the data array. Otherwise, the routine is implementation dependent.
250: - You MUST call VecRestoreArray() when you no longer need access to
251: the array.
252: - Note that the Fortran interface to VecGetArray() differs from the
253: C version. See the users manual for details.
254: */
255: VecGetArray(u,&u_localptr);
257: /*
258: We initialize the solution array by simply writing the solution
259: directly into the array locations. Alternatively, we could use
260: VecSetValues() or VecSetValuesLocal().
261: */
262: for (i=mybase; i<myend; i++) {
263: x = h*(PetscReal)i; /* current location in global grid */
264: u_localptr[i-mybase] = 1.0 + x*x;
265: }
267: /*
268: Restore vector
269: */
270: VecRestoreArray(u,&u_localptr);
272: /*
273: Print debugging information if desired
274: */
275: if (appctx->debug) {
276: PetscPrintf(appctx->comm,"initial guess vector\n");
277: VecView(u,PETSC_VIEWER_STDOUT_WORLD);
278: }
280: return 0;
281: }
282: /* --------------------------------------------------------------------- */
285: /*
286: ExactSolution - Computes the exact solution at a given time.
288: Input Parameters:
289: t - current time
290: solution - vector in which exact solution will be computed
291: appctx - user-defined application context
293: Output Parameter:
294: solution - vector with the newly computed exact solution
295: */
296: PetscErrorCode ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
297: {
298: PetscScalar *s_localptr,h = appctx->h,x;
299: PetscInt i,mybase,myend;
302: /*
303: Determine starting and ending points of each processor's
304: range of grid values
305: */
306: VecGetOwnershipRange(solution,&mybase,&myend);
308: /*
309: Get a pointer to vector data.
310: */
311: VecGetArray(solution,&s_localptr);
313: /*
314: Simply write the solution directly into the array locations.
315: Alternatively, we could use VecSetValues() or VecSetValuesLocal().
316: */
317: for (i=mybase; i<myend; i++) {
318: x = h*(PetscReal)i;
319: s_localptr[i-mybase] = (t + 1.0)*(1.0 + x*x);
320: }
322: /*
323: Restore vector
324: */
325: VecRestoreArray(solution,&s_localptr);
326: return 0;
327: }
328: /* --------------------------------------------------------------------- */
331: /*
332: Monitor - User-provided routine to monitor the solution computed at
333: each timestep. This example plots the solution and computes the
334: error in two different norms.
336: Input Parameters:
337: ts - the timestep context
338: step - the count of the current step (with 0 meaning the
339: initial condition)
340: time - the current time
341: u - the solution at this timestep
342: ctx - the user-provided context for this monitoring routine.
343: In this case we use the application context which contains
344: information about the problem size, workspace and the exact
345: solution.
346: */
347: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec u,void *ctx)
348: {
349: AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */
351: PetscReal en2,en2s,enmax;
352: PetscDraw draw;
354: /*
355: We use the default X windows viewer
356: PETSC_VIEWER_DRAW_(appctx->comm)
357: that is associated with the current communicator. This saves
358: the effort of calling PetscViewerDrawOpen() to create the window.
359: Note that if we wished to plot several items in separate windows we
360: would create each viewer with PetscViewerDrawOpen() and store them in
361: the application context, appctx.
363: PetscReal buffering makes graphics look better.
364: */
365: PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(appctx->comm),0,&draw);
366: PetscDrawSetDoubleBuffer(draw);
367: VecView(u,PETSC_VIEWER_DRAW_(appctx->comm));
369: /*
370: Compute the exact solution at this timestep
371: */
372: ExactSolution(time,appctx->solution,appctx);
374: /*
375: Print debugging information if desired
376: */
377: if (appctx->debug) {
378: PetscPrintf(appctx->comm,"Computed solution vector\n");
379: VecView(u,PETSC_VIEWER_STDOUT_WORLD);
380: PetscPrintf(appctx->comm,"Exact solution vector\n");
381: VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
382: }
384: /*
385: Compute the 2-norm and max-norm of the error
386: */
387: VecAXPY(appctx->solution,-1.0,u);
388: VecNorm(appctx->solution,NORM_2,&en2);
389: en2s = sqrt(appctx->h)*en2; /* scale the 2-norm by the grid spacing */
390: VecNorm(appctx->solution,NORM_MAX,&enmax);
392: /*
393: PetscPrintf() causes only the first processor in this
394: communicator to print the timestep information.
395: */
396: PetscPrintf(appctx->comm,"Timestep %D: time = %G,2-norm error = %G, max norm error = %G\n",
397: step,time,en2s,enmax);
399: /*
400: Print debugging information if desired
401: */
402: if (appctx->debug) {
403: PetscPrintf(appctx->comm,"Error vector\n");
404: VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
405: }
406: return 0;
407: }
408: /* --------------------------------------------------------------------- */
411: /*
412: RHSFunction - User-provided routine that evalues the right-hand-side
413: function of the ODE. This routine is set in the main program by
414: calling TSSetRHSFunction(). We compute:
415: global_out = F(global_in)
417: Input Parameters:
418: ts - timesteping context
419: t - current time
420: global_in - vector containing the current iterate
421: ctx - (optional) user-provided context for function evaluation.
422: In this case we use the appctx defined above.
424: Output Parameter:
425: global_out - vector containing the newly evaluated function
426: */
427: PetscErrorCode RHSFunction(TS ts,PetscReal t,Vec global_in,Vec global_out,void *ctx)
428: {
429: AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */
430: DA da = appctx->da; /* distributed array */
431: Vec local_in = appctx->u_local; /* local ghosted input vector */
432: Vec localwork = appctx->localwork; /* local ghosted work vector */
434: PetscInt i,localsize;
435: PetscMPIInt rank,size;
436: PetscScalar *copyptr,*localptr,sc;
438: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
439: Get ready for local function computations
440: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
441: /*
442: Scatter ghost points to local vector, using the 2-step process
443: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
444: By placing code between these two statements, computations can be
445: done while messages are in transition.
446: */
447: DAGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
448: DAGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);
450: /*
451: Access directly the values in our local INPUT work array
452: */
453: VecGetArray(local_in,&localptr);
455: /*
456: Access directly the values in our local OUTPUT work array
457: */
458: VecGetArray(localwork,©ptr);
460: sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));
462: /*
463: Evaluate our function on the nodes owned by this processor
464: */
465: VecGetLocalSize(local_in,&localsize);
467: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
468: Compute entries for the locally owned part
469: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
471: /*
472: Handle boundary conditions: This is done by using the boundary condition
473: u(t,boundary) = g(t,boundary)
474: for some function g. Now take the derivative with respect to t to obtain
475: u_{t}(t,boundary) = g_{t}(t,boundary)
477: In our case, u(t,0) = t + 1, so that u_{t}(t,0) = 1
478: and u(t,1) = 2t+ 1, so that u_{t}(t,1) = 2
479: */
480: MPI_Comm_rank(appctx->comm,&rank);
481: MPI_Comm_size(appctx->comm,&size);
482: if (!rank) copyptr[0] = 1.0;
483: if (rank == size-1) copyptr[localsize-1] = 2.0;
485: /*
486: Handle the interior nodes where the PDE is replace by finite
487: difference operators.
488: */
489: for (i=1; i<localsize-1; i++) {
490: copyptr[i] = localptr[i] * sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);
491: }
493: /*
494: Restore vectors
495: */
496: VecRestoreArray(local_in,&localptr);
497: VecRestoreArray(localwork,©ptr);
499: /*
500: Insert values from the local OUTPUT vector into the global
501: output vector
502: */
503: DALocalToGlobal(da,localwork,INSERT_VALUES,global_out);
505: /* Print debugging information if desired */
506: if (appctx->debug) {
507: PetscPrintf(appctx->comm,"RHS function vector\n");
508: VecView(global_out,PETSC_VIEWER_STDOUT_WORLD);
509: }
511: return 0;
512: }
513: /* --------------------------------------------------------------------- */
516: /*
517: RHSJacobian - User-provided routine to compute the Jacobian of
518: the nonlinear right-hand-side function of the ODE.
520: Input Parameters:
521: ts - the TS context
522: t - current time
523: global_in - global input vector
524: dummy - optional user-defined context, as set by TSetRHSJacobian()
526: Output Parameters:
527: AA - Jacobian matrix
528: BB - optionally different preconditioning matrix
529: str - flag indicating matrix structure
531: Notes:
532: RHSJacobian computes entries for the locally owned part of the Jacobian.
533: - Currently, all PETSc parallel matrix formats are partitioned by
534: contiguous chunks of rows across the processors.
535: - Each processor needs to insert only elements that it owns
536: locally (but any non-local elements will be sent to the
537: appropriate processor during matrix assembly).
538: - Always specify global row and columns of matrix entries when
539: using MatSetValues().
540: - Here, we set all entries for a particular row at once.
541: - Note that MatSetValues() uses 0-based row and column numbers
542: in Fortran as well as in C.
543: */
544: PetscErrorCode RHSJacobian(TS ts,PetscReal t,Vec global_in,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
545: {
546: Mat A = *AA; /* Jacobian matrix */
547: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
548: Vec local_in = appctx->u_local; /* local ghosted input vector */
549: DA da = appctx->da; /* distributed array */
550: PetscScalar v[3],*localptr,sc;
552: PetscInt i,mstart,mend,mstarts,mends,idx[3],is;
554: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
555: Get ready for local Jacobian computations
556: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
557: /*
558: Scatter ghost points to local vector, using the 2-step process
559: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
560: By placing code between these two statements, computations can be
561: done while messages are in transition.
562: */
563: DAGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
564: DAGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);
566: /*
567: Get pointer to vector data
568: */
569: VecGetArray(local_in,&localptr);
571: /*
572: Get starting and ending locally owned rows of the matrix
573: */
574: MatGetOwnershipRange(A,&mstarts,&mends);
575: mstart = mstarts; mend = mends;
577: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
578: Compute entries for the locally owned part of the Jacobian.
579: - Currently, all PETSc parallel matrix formats are partitioned by
580: contiguous chunks of rows across the processors.
581: - Each processor needs to insert only elements that it owns
582: locally (but any non-local elements will be sent to the
583: appropriate processor during matrix assembly).
584: - Here, we set all entries for a particular row at once.
585: - We can set matrix entries either using either
586: MatSetValuesLocal() or MatSetValues().
587: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
589: /*
590: Set matrix rows corresponding to boundary data
591: */
592: if (mstart == 0) {
593: v[0] = 0.0;
594: MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
595: mstart++;
596: }
597: if (mend == appctx->m) {
598: mend--;
599: v[0] = 0.0;
600: MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
601: }
603: /*
604: Set matrix rows corresponding to interior data. We construct the
605: matrix one row at a time.
606: */
607: sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));
608: for (i=mstart; i<mend; i++) {
609: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
610: is = i - mstart + 1;
611: v[0] = sc*localptr[is];
612: v[1] = sc*(localptr[is+1] + localptr[is-1] - 4.0*localptr[is]);
613: v[2] = sc*localptr[is];
614: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
615: }
617: /*
618: Restore vector
619: */
620: VecRestoreArray(local_in,&localptr);
622: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
623: Complete the matrix assembly process and set some options
624: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
625: /*
626: Assemble matrix, using the 2-step process:
627: MatAssemblyBegin(), MatAssemblyEnd()
628: Computations can be done while messages are in transition
629: by placing code between these two statements.
630: */
631: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
632: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
634: /*
635: Set flag to indicate that the Jacobian matrix retains an identical
636: nonzero structure throughout all timestepping iterations (although the
637: values of the entries change). Thus, we can save some work in setting
638: up the preconditioner (e.g., no need to redo symbolic factorization for
639: ILU/ICC preconditioners).
640: - If the nonzero structure of the matrix is different during
641: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
642: must be used instead. If you are unsure whether the matrix
643: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
644: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
645: believes your assertion and does not check the structure
646: of the matrix. If you erroneously claim that the structure
647: is the same when it actually is not, the new preconditioner
648: will not function correctly. Thus, use this optimization
649: feature with caution!
650: */
651: *str = SAME_NONZERO_PATTERN;
653: /*
654: Set and option to indicate that we will never add a new nonzero location
655: to the matrix. If we do, it will generate an error.
656: */
657: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR);
659: return 0;
660: }