Actual source code: ex3.c
petsc-3.4.2 2013-07-02
2: static char help[] ="Model Equations for Advection-Diffusion\n";
4: /*
5: Page 9, Section 1.2 Model Equations for Advection-Diffusion
7: u_t = a u_x + d u_xx
9: The initial conditions used here different then in the book.
11: */
13: /*
14: Helpful runtime linear solver options:
15: -pc_type mg -da_refine 2 -snes_monitor -ksp_monitor -ts_view (geometric multigrid with three levels)
17: */
19: /*
20: Include "petscts.h" so that we can use TS solvers. Note that this file
21: automatically includes:
22: petscsys.h - base PETSc routines petscvec.h - vectors
23: petscmat.h - matrices
24: petscis.h - index sets petscksp.h - Krylov subspace methods
25: petscviewer.h - viewers petscpc.h - preconditioners
26: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
27: */
29: #include <petscts.h>
30: #include <petscdmda.h>
32: /*
33: User-defined application context - contains data needed by the
34: application-provided call-back routines.
35: */
36: typedef struct {
37: PetscScalar a,d; /* advection and diffusion strength */
38: PetscBool upwind;
39: } AppCtx;
41: /*
42: User-defined routines
43: */
44: extern PetscErrorCode InitialConditions(TS,Vec,AppCtx*);
45: extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Vec,Mat*,Mat*,MatStructure*,void*);
46: extern PetscErrorCode Solution(TS,PetscReal,Vec,AppCtx*);
50: int main(int argc,char **argv)
51: {
52: AppCtx appctx; /* user-defined application context */
53: TS ts; /* timestepping context */
54: Vec U; /* approximate solution vector */
56: PetscReal dt;
57: DM da;
58: PetscInt M;
60: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
61: Initialize program and set problem parameters
62: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
64: PetscInitialize(&argc,&argv,(char*)0,help);
65: appctx.a = 1.0;
66: appctx.d = 0.0;
67: PetscOptionsGetScalar(NULL,"-a",&appctx.a,NULL);
68: PetscOptionsGetScalar(NULL,"-d",&appctx.d,NULL);
69: appctx.upwind = PETSC_TRUE;
70: PetscOptionsGetBool(NULL,"-upwind",&appctx.upwind,NULL);
72: DMDACreate1d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC, -60, 1, 1,NULL,&da);
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Create vector data structures
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: /*
78: Create vector data structures for approximate and exact solutions
79: */
80: DMCreateGlobalVector(da,&U);
82: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
83: Create timestepping solver context
84: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
86: TSCreate(PETSC_COMM_WORLD,&ts);
87: TSSetDM(ts,da);
89: /*
90: For linear problems with a time-dependent f(U,t) in the equation
91: u_t = f(u,t), the user provides the discretized right-hand-side
92: as a time-dependent matrix.
93: */
94: TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);
95: TSSetRHSJacobian(ts,NULL,NULL,RHSMatrixHeat,&appctx);
96: TSSetSolutionFunction(ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void*))Solution,&appctx);
98: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
99: Customize timestepping solver:
100: - Set timestepping duration info
101: Then set runtime options, which can override these defaults.
102: For example,
103: -ts_max_steps <maxsteps> -ts_final_time <maxtime>
104: to override the defaults set by TSSetDuration().
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
108: dt = .48/(M*M);
109: TSSetInitialTimeStep(ts,0.0,dt);
110: TSSetDuration(ts,1000,100.0);
111: TSSetType(ts,TSARKIMEX);
112: TSSetFromOptions(ts);
114: /*
115: Evaluate initial conditions
116: */
117: InitialConditions(ts,U,&appctx);
119: /*
120: Run the timestepping solver
121: */
122: TSSolve(ts,U);
125: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126: Free work space. All PETSc objects should be destroyed when they
127: are no longer needed.
128: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
130: TSDestroy(&ts);
131: VecDestroy(&U);
132: DMDestroy(&da);
134: /*
135: Always call PetscFinalize() before exiting a program. This routine
136: - finalizes the PETSc libraries as well as MPI
137: - provides summary and diagnostic information if certain runtime
138: options are chosen (e.g., -log_summary).
139: */
140: PetscFinalize();
141: return 0;
142: }
143: /* --------------------------------------------------------------------- */
146: /*
147: InitialConditions - Computes the solution at the initial time.
149: Input Parameter:
150: u - uninitialized solution vector (global)
151: appctx - user-defined application context
153: Output Parameter:
154: u - vector with solution at initial time (global)
155: */
156: PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx)
157: {
158: PetscScalar *u,h;
160: PetscInt i,mstart,mend,xm,M;
161: DM da;
163: TSGetDM(ts,&da);
164: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
165: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
166: h = 1.0/M;
167: mend = mstart + xm;
168: /*
169: Get a pointer to vector data.
170: - For default PETSc vectors, VecGetArray() returns a pointer to
171: the data array. Otherwise, the routine is implementation dependent.
172: - You MUST call VecRestoreArray() when you no longer need access to
173: the array.
174: - Note that the Fortran interface to VecGetArray() differs from the
175: C version. See the users manual for details.
176: */
177: DMDAVecGetArray(da,U,&u);
179: /*
180: We initialize the solution array by simply writing the solution
181: directly into the array locations. Alternatively, we could use
182: VecSetValues() or VecSetValuesLocal().
183: */
184: for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
186: /*
187: Restore vector
188: */
189: DMDAVecRestoreArray(da,U,&u);
190: return 0;
191: }
192: /* --------------------------------------------------------------------- */
195: /*
196: Solution - Computes the exact solution at a given time.
198: Input Parameters:
199: t - current time
200: solution - vector in which exact solution will be computed
201: appctx - user-defined application context
203: Output Parameter:
204: solution - vector with the newly computed exact solution
205: */
206: PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx)
207: {
208: PetscScalar *u,ex1,ex2,sc1,sc2,h;
210: PetscInt i,mstart,mend,xm,M;
211: DM da;
213: TSGetDM(ts,&da);
214: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
215: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
216: h = 1.0/M;
217: mend = mstart + xm;
218: /*
219: Get a pointer to vector data.
220: */
221: DMDAVecGetArray(da,U,&u);
223: /*
224: Simply write the solution directly into the array locations.
225: Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
226: */
227: ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*appctx->d*t);
228: ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*appctx->d*t);
229: sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h;
230: for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(sc1*(PetscReal)i + appctx->a*PETSC_PI*6.*t)*ex1 + 3.*PetscSinScalar(sc2*(PetscReal)i + appctx->a*PETSC_PI*2.*t)*ex2;
232: /*
233: Restore vector
234: */
235: DMDAVecRestoreArray(da,U,&u);
236: return 0;
237: }
239: /* --------------------------------------------------------------------- */
242: /*
243: RHSMatrixHeat - User-provided routine to compute the right-hand-side
244: matrix for the heat equation.
246: Input Parameters:
247: ts - the TS context
248: t - current time
249: global_in - global input vector
250: dummy - optional user-defined context, as set by TSetRHSJacobian()
252: Output Parameters:
253: AA - Jacobian matrix
254: BB - optionally different preconditioning matrix
255: str - flag indicating matrix structure
257: Notes:
258: Recall that MatSetValues() uses 0-based row and column numbers
259: in Fortran as well as in C.
260: */
261: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec U,Mat *AA,Mat *BB,MatStructure *str,void *ctx)
262: {
263: Mat A = *AA; /* Jacobian matrix */
264: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
265: PetscInt mstart, mend;
267: PetscInt i,idx[3],M,xm;
268: PetscScalar v[3],h;
269: DM da;
271: TSGetDM(ts,&da);
272: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
273: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
274: h = 1.0/M;
275: mend = mstart + xm;
276: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: Compute entries for the locally owned part of the matrix
278: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279: /*
280: Set matrix rows corresponding to boundary data
281: */
283: /* diffusion */
284: v[0] = appctx->d/(h*h);
285: v[1] = -2.0*appctx->d/(h*h);
286: v[2] = appctx->d/(h*h);
287: if (!mstart) {
288: idx[0] = M-1; idx[1] = 0; idx[2] = 1;
289: MatSetValues(A,1,&mstart,3,idx,v,INSERT_VALUES);
290: mstart++;
291: }
293: if (mend == M) {
294: mend--;
295: idx[0] = M-2; idx[1] = M-1; idx[2] = 0;
296: MatSetValues(A,1,&mend,3,idx,v,INSERT_VALUES);
297: }
299: /*
300: Set matrix rows corresponding to interior data. We construct the
301: matrix one row at a time.
302: */
303: for (i=mstart; i<mend; i++) {
304: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
305: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
306: }
307: MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY);
308: MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY);
310: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
311: mend = mstart + xm;
312: if (!appctx->upwind) {
313: /* advection -- centered differencing */
314: v[0] = -.5*appctx->a/(h);
315: v[1] = .5*appctx->a/(h);
316: if (!mstart) {
317: idx[0] = M-1; idx[1] = 1;
318: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
319: mstart++;
320: }
322: if (mend == M) {
323: mend--;
324: idx[0] = M-2; idx[1] = 0;
325: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
326: }
328: for (i=mstart; i<mend; i++) {
329: idx[0] = i-1; idx[1] = i+1;
330: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
331: }
332: } else {
333: /* advection -- upwinding */
334: v[0] = -appctx->a/(h);
335: v[1] = appctx->a/(h);
336: if (!mstart) {
337: idx[0] = 0; idx[1] = 1;
338: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
339: mstart++;
340: }
342: if (mend == M) {
343: mend--;
344: idx[0] = M-1; idx[1] = 0;
345: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
346: }
348: for (i=mstart; i<mend; i++) {
349: idx[0] = i; idx[1] = i+1;
350: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
351: }
352: }
355: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
356: Complete the matrix assembly process and set some options
357: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
358: /*
359: Assemble matrix, using the 2-step process:
360: MatAssemblyBegin(), MatAssemblyEnd()
361: Computations can be done while messages are in transition
362: by placing code between these two statements.
363: */
364: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
365: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
367: /*
368: Set flag to indicate that the Jacobian matrix retains an identical
369: nonzero structure throughout all timestepping iterations (although the
370: values of the entries change). Thus, we can save some work in setting
371: up the preconditioner (e.g., no need to redo symbolic factorization for
372: ILU/ICC preconditioners).
373: - If the nonzero structure of the matrix is different during
374: successive linear solves, then the flag DIFFERENT_NONZERO_PATTERN
375: must be used instead. If you are unsure whether the matrix
376: structure has changed or not, use the flag DIFFERENT_NONZERO_PATTERN.
377: - Caution: If you specify SAME_NONZERO_PATTERN, PETSc
378: believes your assertion and does not check the structure
379: of the matrix. If you erroneously claim that the structure
380: is the same when it actually is not, the new preconditioner
381: will not function correctly. Thus, use this optimization
382: feature with caution!
383: */
384: *str = SAME_NONZERO_PATTERN;
386: /*
387: Set and option to indicate that we will never add a new nonzero location
388: to the matrix. If we do, it will generate an error.
389: */
390: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
391: return 0;
392: }