Actual source code: ex3.c
3: static char help[] ="Model Equations for Advection-Diffusion\n";
5: /*
6: Page 9, Section 1.2 Model Equations for Advection-Diffusion
8: u_t = a u_x + d u_xx
10: The initial conditions used here different then in the book.
12: */
14: /*
15: Helpful runtime linear solver options:
16: -pc_type mg -da_refine 2 -snes_monitor -ksp_monitor -ts_view (geometric multigrid with three levels)
18: */
20: /*
21: Include "petscts.h" so that we can use TS solvers. Note that this file
22: automatically includes:
23: petscsys.h - base PETSc routines petscvec.h - vectors
24: petscmat.h - matrices
25: petscis.h - index sets petscksp.h - Krylov subspace methods
26: petscviewer.h - viewers petscpc.h - preconditioners
27: petscksp.h - linear solvers petscsnes.h - nonlinear solvers
28: */
30: #include <petscts.h>
31: #include <petscdm.h>
32: #include <petscdmda.h>
34: /*
35: User-defined application context - contains data needed by the
36: application-provided call-back routines.
37: */
38: typedef struct {
39: PetscScalar a,d; /* advection and diffusion strength */
40: PetscBool upwind;
41: } AppCtx;
43: /*
44: User-defined routines
45: */
46: extern PetscErrorCode InitialConditions(TS,Vec,AppCtx*);
47: extern PetscErrorCode RHSMatrixHeat(TS,PetscReal,Vec,Mat,Mat,void*);
48: 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);if (ierr) return ierr;
65: appctx.a = 1.0;
66: appctx.d = 0.0;
67: PetscOptionsGetScalar(NULL,NULL,"-a",&appctx.a,NULL);
68: PetscOptionsGetScalar(NULL,NULL,"-d",&appctx.d,NULL);
69: appctx.upwind = PETSC_TRUE;
70: PetscOptionsGetBool(NULL,NULL,"-upwind",&appctx.upwind,NULL);
72: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);
73: DMSetFromOptions(da);
74: DMSetUp(da);
75: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
76: Create vector data structures
77: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
79: /*
80: Create vector data structures for approximate and exact solutions
81: */
82: DMCreateGlobalVector(da,&U);
84: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
85: Create timestepping solver context
86: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
88: TSCreate(PETSC_COMM_WORLD,&ts);
89: TSSetDM(ts,da);
91: /*
92: For linear problems with a time-dependent f(U,t) in the equation
93: u_t = f(u,t), the user provides the discretized right-hand-side
94: as a time-dependent matrix.
95: */
96: TSSetRHSFunction(ts,NULL,TSComputeRHSFunctionLinear,&appctx);
97: TSSetRHSJacobian(ts,NULL,NULL,RHSMatrixHeat,&appctx);
98: TSSetSolutionFunction(ts,(PetscErrorCode (*)(TS,PetscReal,Vec,void*))Solution,&appctx);
100: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
101: Customize timestepping solver:
102: - Set timestepping duration info
103: Then set runtime options, which can override these defaults.
104: For example,
105: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
106: to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
107: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
109: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
110: dt = .48/(M*M);
111: TSSetTimeStep(ts,dt);
112: TSSetMaxSteps(ts,1000);
113: TSSetMaxTime(ts,100.0);
114: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
115: TSSetType(ts,TSARKIMEX);
116: TSSetFromOptions(ts);
118: /*
119: Evaluate initial conditions
120: */
121: InitialConditions(ts,U,&appctx);
123: /*
124: Run the timestepping solver
125: */
126: TSSolve(ts,U);
129: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130: Free work space. All PETSc objects should be destroyed when they
131: are no longer needed.
132: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
134: TSDestroy(&ts);
135: VecDestroy(&U);
136: DMDestroy(&da);
138: /*
139: Always call PetscFinalize() before exiting a program. This routine
140: - finalizes the PETSc libraries as well as MPI
141: - provides summary and diagnostic information if certain runtime
142: options are chosen (e.g., -log_view).
143: */
144: PetscFinalize();
145: return ierr;
146: }
147: /* --------------------------------------------------------------------- */
148: /*
149: InitialConditions - Computes the solution at the initial time.
151: Input Parameter:
152: u - uninitialized solution vector (global)
153: appctx - user-defined application context
155: Output Parameter:
156: u - vector with solution at initial time (global)
157: */
158: PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx)
159: {
160: PetscScalar *u,h;
162: PetscInt i,mstart,mend,xm,M;
163: DM da;
165: TSGetDM(ts,&da);
166: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
167: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
168: h = 1.0/M;
169: mend = mstart + xm;
170: /*
171: Get a pointer to vector data.
172: - For default PETSc vectors, VecGetArray() returns a pointer to
173: the data array. Otherwise, the routine is implementation dependent.
174: - You MUST call VecRestoreArray() when you no longer need access to
175: the array.
176: - Note that the Fortran interface to VecGetArray() differs from the
177: C version. See the users manual for details.
178: */
179: DMDAVecGetArray(da,U,&u);
181: /*
182: We initialize the solution array by simply writing the solution
183: directly into the array locations. Alternatively, we could use
184: VecSetValues() or VecSetValuesLocal().
185: */
186: for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
188: /*
189: Restore vector
190: */
191: DMDAVecRestoreArray(da,U,&u);
192: return 0;
193: }
194: /* --------------------------------------------------------------------- */
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: /* --------------------------------------------------------------------- */
240: /*
241: RHSMatrixHeat - User-provided routine to compute the right-hand-side
242: matrix for the heat equation.
244: Input Parameters:
245: ts - the TS context
246: t - current time
247: global_in - global input vector
248: dummy - optional user-defined context, as set by TSetRHSJacobian()
250: Output Parameters:
251: AA - Jacobian matrix
252: BB - optionally different preconditioning matrix
253: str - flag indicating matrix structure
255: Notes:
256: Recall that MatSetValues() uses 0-based row and column numbers
257: in Fortran as well as in C.
258: */
259: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec U,Mat AA,Mat BB,void *ctx)
260: {
261: Mat A = AA; /* Jacobian matrix */
262: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
263: PetscInt mstart, mend;
265: PetscInt i,idx[3],M,xm;
266: PetscScalar v[3],h;
267: DM da;
269: TSGetDM(ts,&da);
270: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
271: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
272: h = 1.0/M;
273: mend = mstart + xm;
274: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
275: Compute entries for the locally owned part of the matrix
276: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
277: /*
278: Set matrix rows corresponding to boundary data
279: */
281: /* diffusion */
282: v[0] = appctx->d/(h*h);
283: v[1] = -2.0*appctx->d/(h*h);
284: v[2] = appctx->d/(h*h);
285: if (!mstart) {
286: idx[0] = M-1; idx[1] = 0; idx[2] = 1;
287: MatSetValues(A,1,&mstart,3,idx,v,INSERT_VALUES);
288: mstart++;
289: }
291: if (mend == M) {
292: mend--;
293: idx[0] = M-2; idx[1] = M-1; idx[2] = 0;
294: MatSetValues(A,1,&mend,3,idx,v,INSERT_VALUES);
295: }
297: /*
298: Set matrix rows corresponding to interior data. We construct the
299: matrix one row at a time.
300: */
301: for (i=mstart; i<mend; i++) {
302: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
303: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
304: }
305: MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY);
306: MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY);
308: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
309: mend = mstart + xm;
310: if (!appctx->upwind) {
311: /* advection -- centered differencing */
312: v[0] = -.5*appctx->a/(h);
313: v[1] = .5*appctx->a/(h);
314: if (!mstart) {
315: idx[0] = M-1; idx[1] = 1;
316: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
317: mstart++;
318: }
320: if (mend == M) {
321: mend--;
322: idx[0] = M-2; idx[1] = 0;
323: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
324: }
326: for (i=mstart; i<mend; i++) {
327: idx[0] = i-1; idx[1] = i+1;
328: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
329: }
330: } else {
331: /* advection -- upwinding */
332: v[0] = -appctx->a/(h);
333: v[1] = appctx->a/(h);
334: if (!mstart) {
335: idx[0] = 0; idx[1] = 1;
336: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
337: mstart++;
338: }
340: if (mend == M) {
341: mend--;
342: idx[0] = M-1; idx[1] = 0;
343: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
344: }
346: for (i=mstart; i<mend; i++) {
347: idx[0] = i; idx[1] = i+1;
348: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
349: }
350: }
353: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
354: Complete the matrix assembly process and set some options
355: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
356: /*
357: Assemble matrix, using the 2-step process:
358: MatAssemblyBegin(), MatAssemblyEnd()
359: Computations can be done while messages are in transition
360: by placing code between these two statements.
361: */
362: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
363: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
365: /*
366: Set and option to indicate that we will never add a new nonzero location
367: to the matrix. If we do, it will generate an error.
368: */
369: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
370: return 0;
371: }
374: /*TEST
376: test:
377: args: -pc_type mg -da_refine 2 -ts_view -ts_monitor -ts_max_time .3 -mg_levels_ksp_max_it 3
378: requires: double
379: filter: grep -v "total number of"
381: test:
382: suffix: 2
383: 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
384: requires: x
385: output_file: output/ex3_1.out
386: requires: double
387: filter: grep -v "total number of"
389: TEST*/