Actual source code: ex3.c
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 <petscdm.h>
31: #include <petscdmda.h>
33: /*
34: User-defined application context - contains data needed by the
35: application-provided call-back routines.
36: */
37: typedef struct {
38: PetscScalar a,d; /* advection and diffusion strength */
39: PetscBool upwind;
40: } AppCtx;
42: /*
43: User-defined routines
44: */
45: extern PetscErrorCode InitialConditions(TS,Vec,AppCtx*);
46: extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Vec,Mat,Mat,void*);
47: extern PetscErrorCode Solution(TS,PetscReal,Vec,AppCtx*);
49: int main(int argc,char **argv)
50: {
51: AppCtx appctx; /* user-defined application context */
52: TS ts; /* timestepping context */
53: Vec U; /* approximate solution vector */
55: PetscReal dt;
56: DM da;
57: PetscInt M;
59: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
60: Initialize program and set problem parameters
61: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
63: PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
64: appctx.a = 1.0;
65: appctx.d = 0.0;
66: PetscOptionsGetScalar(NULL,NULL,"-a",&appctx.a,NULL);
67: PetscOptionsGetScalar(NULL,NULL,"-d",&appctx.d,NULL);
68: appctx.upwind = PETSC_TRUE;
69: PetscOptionsGetBool(NULL,NULL,"-upwind",&appctx.upwind,NULL);
71: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);
72: DMSetFromOptions(da);
73: DMSetUp(da);
74: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
75: Create vector data structures
76: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
78: /*
79: Create vector data structures for approximate and exact solutions
80: */
81: DMCreateGlobalVector(da,&U);
83: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
84: Create timestepping solver context
85: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
87: TSCreate(PETSC_COMM_WORLD,&ts);
88: TSSetDM(ts,da);
90: /*
91: For linear problems with a time-dependent f(U,t) in the equation
92: u_t = f(u,t), the user provides the discretized right-hand-side
93: as a time-dependent matrix.
94: */
95: TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);
96: TSSetRHSJacobian(ts,NULL,NULL,RHSMatrixHeat,&appctx);
97: TSSetSolutionFunction(ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void*))Solution,&appctx);
99: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
100: Customize timestepping solver:
101: - Set timestepping duration info
102: Then set runtime options, which can override these defaults.
103: For example,
104: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
105: to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
106: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
108: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
109: dt = .48/(M*M);
110: TSSetTimeStep(ts,dt);
111: TSSetMaxSteps(ts,1000);
112: TSSetMaxTime(ts,100.0);
113: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
114: TSSetType(ts,TSARKIMEX);
115: TSSetFromOptions(ts);
117: /*
118: Evaluate initial conditions
119: */
120: InitialConditions(ts,U,&appctx);
122: /*
123: Run the timestepping solver
124: */
125: TSSolve(ts,U);
127: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
128: Free work space. All PETSc objects should be destroyed when they
129: are no longer needed.
130: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132: TSDestroy(&ts);
133: VecDestroy(&U);
134: DMDestroy(&da);
136: /*
137: Always call PetscFinalize() before exiting a program. This routine
138: - finalizes the PETSc libraries as well as MPI
139: - provides summary and diagnostic information if certain runtime
140: options are chosen (e.g., -log_view).
141: */
142: PetscFinalize();
143: return ierr;
144: }
145: /* --------------------------------------------------------------------- */
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: /* --------------------------------------------------------------------- */
193: /*
194: Solution - Computes the exact solution at a given time.
196: Input Parameters:
197: t - current time
198: solution - vector in which exact solution will be computed
199: appctx - user-defined application context
201: Output Parameter:
202: solution - vector with the newly computed exact solution
203: */
204: PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx)
205: {
206: PetscScalar *u,ex1,ex2,sc1,sc2,h;
208: PetscInt i,mstart,mend,xm,M;
209: DM da;
211: TSGetDM(ts,&da);
212: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
213: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
214: h = 1.0/M;
215: mend = mstart + xm;
216: /*
217: Get a pointer to vector data.
218: */
219: DMDAVecGetArray(da,U,&u);
221: /*
222: Simply write the solution directly into the array locations.
223: Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
224: */
225: ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*appctx->d*t);
226: ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*appctx->d*t);
227: sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h;
228: 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;
230: /*
231: Restore vector
232: */
233: DMDAVecRestoreArray(da,U,&u);
234: return 0;
235: }
237: /* --------------------------------------------------------------------- */
238: /*
239: RHSMatrixHeat - User-provided routine to compute the right-hand-side
240: matrix for the heat equation.
242: Input Parameters:
243: ts - the TS context
244: t - current time
245: global_in - global input vector
246: dummy - optional user-defined context, as set by TSetRHSJacobian()
248: Output Parameters:
249: AA - Jacobian matrix
250: BB - optionally different preconditioning matrix
251: str - flag indicating matrix structure
253: Notes:
254: Recall that MatSetValues() uses 0-based row and column numbers
255: in Fortran as well as in C.
256: */
257: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec U,Mat AA,Mat BB,void *ctx)
258: {
259: Mat A = AA; /* Jacobian matrix */
260: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
261: PetscInt mstart, mend;
263: PetscInt i,idx[3],M,xm;
264: PetscScalar v[3],h;
265: DM da;
267: TSGetDM(ts,&da);
268: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
269: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
270: h = 1.0/M;
271: mend = mstart + xm;
272: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
273: Compute entries for the locally owned part of the matrix
274: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
275: /*
276: Set matrix rows corresponding to boundary data
277: */
279: /* diffusion */
280: v[0] = appctx->d/(h*h);
281: v[1] = -2.0*appctx->d/(h*h);
282: v[2] = appctx->d/(h*h);
283: if (!mstart) {
284: idx[0] = M-1; idx[1] = 0; idx[2] = 1;
285: MatSetValues(A,1,&mstart,3,idx,v,INSERT_VALUES);
286: mstart++;
287: }
289: if (mend == M) {
290: mend--;
291: idx[0] = M-2; idx[1] = M-1; idx[2] = 0;
292: MatSetValues(A,1,&mend,3,idx,v,INSERT_VALUES);
293: }
295: /*
296: Set matrix rows corresponding to interior data. We construct the
297: matrix one row at a time.
298: */
299: for (i=mstart; i<mend; i++) {
300: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
301: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
302: }
303: MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY);
304: MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY);
306: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
307: mend = mstart + xm;
308: if (!appctx->upwind) {
309: /* advection -- centered differencing */
310: v[0] = -.5*appctx->a/(h);
311: v[1] = .5*appctx->a/(h);
312: if (!mstart) {
313: idx[0] = M-1; idx[1] = 1;
314: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
315: mstart++;
316: }
318: if (mend == M) {
319: mend--;
320: idx[0] = M-2; idx[1] = 0;
321: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
322: }
324: for (i=mstart; i<mend; i++) {
325: idx[0] = i-1; idx[1] = i+1;
326: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
327: }
328: } else {
329: /* advection -- upwinding */
330: v[0] = -appctx->a/(h);
331: v[1] = appctx->a/(h);
332: if (!mstart) {
333: idx[0] = 0; idx[1] = 1;
334: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
335: mstart++;
336: }
338: if (mend == M) {
339: mend--;
340: idx[0] = M-1; idx[1] = 0;
341: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
342: }
344: for (i=mstart; i<mend; i++) {
345: idx[0] = i; idx[1] = i+1;
346: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
347: }
348: }
350: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
351: Complete the matrix assembly process and set some options
352: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
353: /*
354: Assemble matrix, using the 2-step process:
355: MatAssemblyBegin(), MatAssemblyEnd()
356: Computations can be done while messages are in transition
357: by placing code between these two statements.
358: */
359: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
360: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
362: /*
363: Set and option to indicate that we will never add a new nonzero location
364: to the matrix. If we do, it will generate an error.
365: */
366: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
367: return 0;
368: }
370: /*TEST
372: test:
373: args: -pc_type mg -da_refine 2 -ts_view -ts_monitor -ts_max_time .3 -mg_levels_ksp_max_it 3
374: requires: double
375: filter: grep -v "total number of"
377: test:
378: suffix: 2
379: args: -pc_type mg -da_refine 2 -ts_view -ts_monitor_draw_solution -ts_monitor -ts_max_time .3 -mg_levels_ksp_max_it 3
380: requires: x
381: output_file: output/ex3_1.out
382: requires: double
383: filter: grep -v "total number of"
385: TEST*/