Actual source code: ex2.c
1: /*$Id: ex2.c,v 1.41 2001/08/10 03:34:17 bsmith Exp $*/
2: static char help[] ="Solves a time-dependent nonlinear PDE. Uses implicitn
3: timestepping. Runtime options include:n
4: -M <xg>, where <xg> = number of grid pointsn
5: -debug : Activate debugging printoutsn
6: -nox : Deactivate x-window graphicsnn";
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: int 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: */
65: extern int InitialConditions(Vec,AppCtx*);
66: extern int RHSFunction(TS,PetscReal,Vec,Vec,void*);
67: extern int RHSJacobian(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
68: extern int Monitor(TS,int,PetscReal,Vec,void*);
69: extern int ExactSolution(PetscReal,Vec,AppCtx*);
71: /*
72: Utility routine for finite difference Jacobian approximation
73: */
74: extern int RHSJacobianFD(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
76: #undef __FUNCT__
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: int time_steps_max = 1000; /* default max timesteps */
85: int ierr,steps;
86: PetscReal ftime; /* final time */
87: PetscReal dt;
88: PetscReal time_total_max = 100.0; /* default max total time */
89: PetscTruth flg;
91: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
92: Initialize program and set problem parameters
93: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
94:
95: PetscInitialize(&argc,&argv,(char*)0,help);
97: appctx.comm = PETSC_COMM_WORLD;
98: appctx.m = 60;
99: PetscOptionsGetInt(PETSC_NULL,"-M",&appctx.m,PETSC_NULL);
100: PetscOptionsHasName(PETSC_NULL,"-debug",&appctx.debug);
101: appctx.h = 1.0/(appctx.m-1.0);
103: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
104: Create vector data structures
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: /*
108: Create distributed array (DA) to manage parallel grid and vectors
109: and to set up the ghost point communication pattern. There are M
110: total grid values spread equally among all the processors.
111: */
112: DACreate1d(PETSC_COMM_WORLD,DA_NONPERIODIC,appctx.m,1,1,PETSC_NULL,
113: &appctx.da);
115: /*
116: Extract global and local vectors from DA; we use these to store the
117: approximate solution. Then duplicate these for remaining vectors that
118: have the same types.
119: */
120: DACreateGlobalVector(appctx.da,&u);
121: DACreateLocalVector(appctx.da,&appctx.u_local);
123: /*
124: Create local work vector for use in evaluating right-hand-side function;
125: create global work vector for storing exact solution.
126: */
127: VecDuplicate(appctx.u_local,&appctx.localwork);
128: VecDuplicate(u,&appctx.solution);
130: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
131: Create timestepping solver context; set callback routine for
132: right-hand-side function evaluation.
133: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
135: TSCreate(PETSC_COMM_WORLD,&ts);
136: TSSetProblemType(ts,TS_NONLINEAR);
137: TSSetRHSFunction(ts,RHSFunction,&appctx);
139: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
140: Set optional user-defined monitoring routine
141: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
143: TSSetMonitor(ts,Monitor,&appctx,PETSC_NULL);
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: For nonlinear problems, the user can provide a Jacobian evaluation
147: routine (or use a finite differencing approximation).
149: Create matrix data structure; set Jacobian evaluation routine.
150: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152: MatCreate(PETSC_COMM_WORLD,PETSC_DECIDE,PETSC_DECIDE,appctx.m,appctx.m,&A);
153: MatSetFromOptions(A);
154: PetscOptionsHasName(PETSC_NULL,"-fdjac",&flg);
155: if (flg) {
156: TSSetRHSJacobian(ts,A,A,RHSJacobianFD,&appctx);
157: } else {
158: TSSetRHSJacobian(ts,A,A,RHSJacobian,&appctx);
159: }
161: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
162: Set solution vector and initial timestep
163: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
165: dt = appctx.h/2.0;
166: TSSetInitialTimeStep(ts,0.0,dt);
167: TSSetSolution(ts,u);
169: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
170: Customize timestepping solver:
171: - Set the solution method to be the Backward Euler method.
172: - Set timestepping duration info
173: Then set runtime options, which can override these defaults.
174: For example,
175: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
176: to override the defaults set by TSSetDuration().
177: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179: TSSetType(ts,TS_BEULER);
180: TSSetDuration(ts,time_steps_max,time_total_max);
181: TSSetFromOptions(ts);
183: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
184: Solve the problem
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187: /*
188: Evaluate initial conditions
189: */
190: InitialConditions(u,&appctx);
192: /*
193: Run the timestepping solver
194: */
195: TSStep(ts,&steps,&ftime);
197: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
198: Free work space. All PETSc objects should be destroyed when they
199: are no longer needed.
200: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
202: TSDestroy(ts);
203: VecDestroy(u);
204: MatDestroy(A);
205: DADestroy(appctx.da);
206: VecDestroy(appctx.localwork);
207: VecDestroy(appctx.solution);
208: VecDestroy(appctx.u_local);
210: /*
211: Always call PetscFinalize() before exiting a program. This routine
212: - finalizes the PETSc libraries as well as MPI
213: - provides summary and diagnostic information if certain runtime
214: options are chosen (e.g., -log_summary).
215: */
216: PetscFinalize();
217: return 0;
218: }
219: /* --------------------------------------------------------------------- */
220: #undef __FUNCT__
222: /*
223: InitialConditions - Computes the solution at the initial time.
225: Input Parameters:
226: u - uninitialized solution vector (global)
227: appctx - user-defined application context
229: Output Parameter:
230: u - vector with solution at initial time (global)
231: */
232: int InitialConditions(Vec u,AppCtx *appctx)
233: {
234: PetscScalar *u_localptr,h = appctx->h,x;
235: int i,mybase,myend,ierr;
237: /*
238: Determine starting point of each processor's range of
239: grid values.
240: */
241: VecGetOwnershipRange(u,&mybase,&myend);
243: /*
244: Get a pointer to vector data.
245: - For default PETSc vectors, VecGetArray() returns a pointer to
246: the data array. Otherwise, the routine is implementation dependent.
247: - You MUST call VecRestoreArray() when you no longer need access to
248: the array.
249: - Note that the Fortran interface to VecGetArray() differs from the
250: C version. See the users manual for details.
251: */
252: VecGetArray(u,&u_localptr);
254: /*
255: We initialize the solution array by simply writing the solution
256: directly into the array locations. Alternatively, we could use
257: VecSetValues() or VecSetValuesLocal().
258: */
259: for (i=mybase; i<myend; i++) {
260: x = h*(PetscReal)i; /* current location in global grid */
261: u_localptr[i-mybase] = 1.0 + x*x;
262: }
264: /*
265: Restore vector
266: */
267: VecRestoreArray(u,&u_localptr);
269: /*
270: Print debugging information if desired
271: */
272: if (appctx->debug) {
273: PetscPrintf(appctx->comm,"initial guess vectorn");
274: VecView(u,PETSC_VIEWER_STDOUT_WORLD);
275: }
277: return 0;
278: }
279: /* --------------------------------------------------------------------- */
280: #undef __FUNCT__
282: /*
283: ExactSolution - Computes the exact solution at a given time.
285: Input Parameters:
286: t - current time
287: solution - vector in which exact solution will be computed
288: appctx - user-defined application context
290: Output Parameter:
291: solution - vector with the newly computed exact solution
292: */
293: int ExactSolution(PetscReal t,Vec solution,AppCtx *appctx)
294: {
295: PetscScalar *s_localptr,h = appctx->h,x;
296: int i,mybase,myend,ierr;
298: /*
299: Determine starting and ending points of each processor's
300: range of grid values
301: */
302: VecGetOwnershipRange(solution,&mybase,&myend);
304: /*
305: Get a pointer to vector data.
306: */
307: VecGetArray(solution,&s_localptr);
309: /*
310: Simply write the solution directly into the array locations.
311: Alternatively, we could use VecSetValues() or VecSetValuesLocal().
312: */
313: for (i=mybase; i<myend; i++) {
314: x = h*(PetscReal)i;
315: s_localptr[i-mybase] = (t + 1.0)*(1.0 + x*x);
316: }
318: /*
319: Restore vector
320: */
321: VecRestoreArray(solution,&s_localptr);
322: return 0;
323: }
324: /* --------------------------------------------------------------------- */
325: #undef __FUNCT__
327: /*
328: Monitor - User-provided routine to monitor the solution computed at
329: each timestep. This example plots the solution and computes the
330: error in two different norms.
332: Input Parameters:
333: ts - the timestep context
334: step - the count of the current step (with 0 meaning the
335: initial condition)
336: time - the current time
337: u - the solution at this timestep
338: ctx - the user-provided context for this monitoring routine.
339: In this case we use the application context which contains
340: information about the problem size, workspace and the exact
341: solution.
342: */
343: int Monitor(TS ts,int step,PetscReal time,Vec u,void *ctx)
344: {
345: AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */
346: int ierr;
347: PetscReal en2,en2s,enmax;
348: PetscScalar mone = -1.0;
349: PetscDraw draw;
351: /*
352: We use the default X windows viewer
353: PETSC_VIEWER_DRAW_(appctx->comm)
354: that is associated with the current communicator. This saves
355: the effort of calling PetscViewerDrawOpen() to create the window.
356: Note that if we wished to plot several items in separate windows we
357: would create each viewer with PetscViewerDrawOpen() and store them in
358: the application context, appctx.
360: PetscReal buffering makes graphics look better.
361: */
362: PetscViewerDrawGetDraw(PETSC_VIEWER_DRAW_(appctx->comm),0,&draw);
363: PetscDrawSetDoubleBuffer(draw);
364: VecView(u,PETSC_VIEWER_DRAW_(appctx->comm));
366: /*
367: Compute the exact solution at this timestep
368: */
369: ExactSolution(time,appctx->solution,appctx);
371: /*
372: Print debugging information if desired
373: */
374: if (appctx->debug) {
375: PetscPrintf(appctx->comm,"Computed solution vectorn");
376: VecView(u,PETSC_VIEWER_STDOUT_WORLD);
377: PetscPrintf(appctx->comm,"Exact solution vectorn");
378: VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
379: }
381: /*
382: Compute the 2-norm and max-norm of the error
383: */
384: VecAXPY(&mone,u,appctx->solution);
385: VecNorm(appctx->solution,NORM_2,&en2);
386: en2s = sqrt(appctx->h)*en2; /* scale the 2-norm by the grid spacing */
387: VecNorm(appctx->solution,NORM_MAX,&enmax);
389: /*
390: PetscPrintf() causes only the first processor in this
391: communicator to print the timestep information.
392: */
393: PetscPrintf(appctx->comm,"Timestep %d: time = %g,2-norm error = %g, max norm error = %gn",
394: step,time,en2s,enmax);
396: /*
397: Print debugging information if desired
398: */
399: if (appctx->debug) {
400: PetscPrintf(appctx->comm,"Error vectorn");
401: VecView(appctx->solution,PETSC_VIEWER_STDOUT_WORLD);
402: }
403: return 0;
404: }
405: /* --------------------------------------------------------------------- */
406: #undef __FUNCT__
408: /*
409: RHSFunction - User-provided routine that evalues the right-hand-side
410: function of the ODE. This routine is set in the main program by
411: calling TSSetRHSFunction(). We compute:
412: global_out = F(global_in)
414: Input Parameters:
415: ts - timesteping context
416: t - current time
417: global_in - vector containing the current iterate
418: ctx - (optional) user-provided context for function evaluation.
419: In this case we use the appctx defined above.
421: Output Parameter:
422: global_out - vector containing the newly evaluated function
423: */
424: int RHSFunction(TS ts,PetscReal t,Vec global_in,Vec global_out,void *ctx)
425: {
426: AppCtx *appctx = (AppCtx*) ctx; /* user-defined application context */
427: DA da = appctx->da; /* distributed array */
428: Vec local_in = appctx->u_local; /* local ghosted input vector */
429: Vec localwork = appctx->localwork; /* local ghosted work vector */
430: int ierr,i,localsize,rank,size;
431: PetscScalar *copyptr,*localptr,sc;
433: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
434: Get ready for local function computations
435: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
436: /*
437: Scatter ghost points to local vector, using the 2-step process
438: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
439: By placing code between these two statements, computations can be
440: done while messages are in transition.
441: */
442: DAGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
443: DAGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);
445: /*
446: Access directly the values in our local INPUT work array
447: */
448: VecGetArray(local_in,&localptr);
450: /*
451: Access directly the values in our local OUTPUT work array
452: */
453: VecGetArray(localwork,©ptr);
455: sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));
457: /*
458: Evaluate our function on the nodes owned by this processor
459: */
460: VecGetLocalSize(local_in,&localsize);
462: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
463: Compute entries for the locally owned part
464: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
466: /*
467: Handle boundary conditions: This is done by using the boundary condition
468: u(t,boundary) = g(t,boundary)
469: for some function g. Now take the derivative with respect to t to obtain
470: u_{t}(t,boundary) = g_{t}(t,boundary)
472: In our case, u(t,0) = t + 1, so that u_{t}(t,0) = 1
473: and u(t,1) = 2t+ 1, so that u_{t}(t,1) = 2
474: */
475: MPI_Comm_rank(appctx->comm,&rank);
476: MPI_Comm_size(appctx->comm,&size);
477: if (!rank) copyptr[0] = 1.0;
478: if (rank == size-1) copyptr[localsize-1] = 2.0;
480: /*
481: Handle the interior nodes where the PDE is replace by finite
482: difference operators.
483: */
484: for (i=1; i<localsize-1; i++) {
485: copyptr[i] = localptr[i] * sc * (localptr[i+1] + localptr[i-1] - 2.0*localptr[i]);
486: }
488: /*
489: Restore vectors
490: */
491: VecRestoreArray(local_in,&localptr);
492: VecRestoreArray(localwork,©ptr);
494: /*
495: Insert values from the local OUTPUT vector into the global
496: output vector
497: */
498: DALocalToGlobal(da,localwork,INSERT_VALUES,global_out);
500: /* Print debugging information if desired */
501: if (appctx->debug) {
502: PetscPrintf(appctx->comm,"RHS function vectorn");
503: VecView(global_out,PETSC_VIEWER_STDOUT_WORLD);
504: }
506: return 0;
507: }
508: /* --------------------------------------------------------------------- */
509: #undef __FUNCT__
511: /*
512: RHSJacobian - User-provided routine to compute the Jacobian of
513: the nonlinear right-hand-side function of the ODE.
515: Input Parameters:
516: ts - the TS context
517: t - current time
518: global_in - global input vector
519: dummy - optional user-defined context, as set by TSetRHSJacobian()
521: Output Parameters:
522: AA - Jacobian matrix
523: BB - optionally different preconditioning matrix
524: str - flag indicating matrix structure
526: Notes:
527: RHSJacobian computes entries for the locally owned part of the Jacobian.
528: - Currently, all PETSc parallel matrix formats are partitioned by
529: contiguous chunks of rows across the processors.
530: - Each processor needs to insert only elements that it owns
531: locally (but any non-local elements will be sent to the
532: appropriate processor during matrix assembly).
533: - Always specify global row and columns of matrix entries when
534: using MatSetValues().
535: - Here, we set all entries for a particular row at once.
536: - Note that MatSetValues() uses 0-based row and column numbers
537: in Fortran as well as in C.
538: */
539: int RHSJacobian(TS ts,PetscReal t,Vec global_in,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
540: {
541: Mat A = *AA; /* Jacobian matrix */
542: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
543: Vec local_in = appctx->u_local; /* local ghosted input vector */
544: DA da = appctx->da; /* distributed array */
545: PetscScalar v[3],*localptr,sc;
546: int ierr,i,mstart,mend,mstarts,mends,idx[3],is;
548: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
549: Get ready for local Jacobian computations
550: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
551: /*
552: Scatter ghost points to local vector, using the 2-step process
553: DAGlobalToLocalBegin(), DAGlobalToLocalEnd().
554: By placing code between these two statements, computations can be
555: done while messages are in transition.
556: */
557: DAGlobalToLocalBegin(da,global_in,INSERT_VALUES,local_in);
558: DAGlobalToLocalEnd(da,global_in,INSERT_VALUES,local_in);
560: /*
561: Get pointer to vector data
562: */
563: VecGetArray(local_in,&localptr);
565: /*
566: Get starting and ending locally owned rows of the matrix
567: */
568: MatGetOwnershipRange(A,&mstarts,&mends);
569: mstart = mstarts; mend = mends;
571: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
572: Compute entries for the locally owned part of the Jacobian.
573: - Currently, all PETSc parallel matrix formats are partitioned by
574: contiguous chunks of rows across the processors.
575: - Each processor needs to insert only elements that it owns
576: locally (but any non-local elements will be sent to the
577: appropriate processor during matrix assembly).
578: - Here, we set all entries for a particular row at once.
579: - We can set matrix entries either using either
580: MatSetValuesLocal() or MatSetValues().
581: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
583: /*
584: Set matrix rows corresponding to boundary data
585: */
586: if (mstart == 0) {
587: v[0] = 0.0;
588: MatSetValues(A,1,&mstart,1,&mstart,v,INSERT_VALUES);
589: mstart++;
590: }
591: if (mend == appctx->m) {
592: mend--;
593: v[0] = 0.0;
594: MatSetValues(A,1,&mend,1,&mend,v,INSERT_VALUES);
595: }
597: /*
598: Set matrix rows corresponding to interior data. We construct the
599: matrix one row at a time.
600: */
601: sc = 1.0/(appctx->h*appctx->h*2.0*(1.0+t)*(1.0+t));
602: for (i=mstart; i<mend; i++) {
603: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
604: is = i - mstart + 1;
605: v[0] = sc*localptr[is];
606: v[1] = sc*(localptr[is+1] + localptr[is-1] - 4.0*localptr[is]);
607: v[2] = sc*localptr[is];
608: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
609: }
611: /*
612: Restore vector
613: */
614: VecRestoreArray(local_in,&localptr);
616: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
617: Complete the matrix assembly process and set some options
618: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
619: /*
620: Assemble matrix, using the 2-step process:
621: MatAssemblyBegin(), MatAssemblyEnd()
622: Computations can be done while messages are in transition
623: by placing code between these two statements.
624: */
625: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
626: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
628: /*
629: Set flag to indicate that the Jacobian matrix retains an identical
630: nonzero structure throughout all timestepping iterations (although the
631: values of the entries change). Thus, we can save some work in setting
632: up the preconditioner (e.g., no need to redo symbolic factorization for
633: ILU/ICC preconditioners).
634: - If the nonzero structure of the matrix is different during
635: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
636: must be used instead. If you are unsure whether the matrix
637: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
638: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
639: believes your assertion and does not check the structure
640: of the matrix. If you erroneously claim that the structure
641: is the same when it actually is not, the new preconditioner
642: will not function correctly. Thus, use this optimization
643: feature with caution!
644: */
645: *str = SAME_NONZERO_PATTERN;
647: /*
648: Set and option to indicate that we will never add a new nonzero location
649: to the matrix. If we do, it will generate an error.
650: */
651: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR);
653: return 0;
654: }