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 */
54: PetscReal dt;
55: DM da;
56: PetscInt M;
58: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
59: Initialize program and set problem parameters
60: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
62: PetscInitialize(&argc,&argv,(char*)0,help);
63: appctx.a = 1.0;
64: appctx.d = 0.0;
65: PetscOptionsGetScalar(NULL,NULL,"-a",&appctx.a,NULL);
66: PetscOptionsGetScalar(NULL,NULL,"-d",&appctx.d,NULL);
67: appctx.upwind = PETSC_TRUE;
68: PetscOptionsGetBool(NULL,NULL,"-upwind",&appctx.upwind,NULL);
70: DMDACreate1d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC, 60, 1, 1,NULL,&da);
71: DMSetFromOptions(da);
72: DMSetUp(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_max_time <maxtime>
104: to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
105: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
107: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
108: dt = .48/(M*M);
109: TSSetTimeStep(ts,dt);
110: TSSetMaxSteps(ts,1000);
111: TSSetMaxTime(ts,100.0);
112: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);
113: TSSetType(ts,TSARKIMEX);
114: TSSetFromOptions(ts);
116: /*
117: Evaluate initial conditions
118: */
119: InitialConditions(ts,U,&appctx);
121: /*
122: Run the timestepping solver
123: */
124: TSSolve(ts,U);
126: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127: Free work space. All PETSc objects should be destroyed when they
128: are no longer needed.
129: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
131: TSDestroy(&ts);
132: VecDestroy(&U);
133: DMDestroy(&da);
135: /*
136: Always call PetscFinalize() before exiting a program. This routine
137: - finalizes the PETSc libraries as well as MPI
138: - provides summary and diagnostic information if certain runtime
139: options are chosen (e.g., -log_view).
140: */
141: PetscFinalize();
142: return 0;
143: }
144: /* --------------------------------------------------------------------- */
145: /*
146: InitialConditions - Computes the solution at the initial time.
148: Input Parameter:
149: u - uninitialized solution vector (global)
150: appctx - user-defined application context
152: Output Parameter:
153: u - vector with solution at initial time (global)
154: */
155: PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx)
156: {
157: PetscScalar *u,h;
158: PetscInt i,mstart,mend,xm,M;
159: DM da;
161: TSGetDM(ts,&da);
162: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
163: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
164: h = 1.0/M;
165: mend = mstart + xm;
166: /*
167: Get a pointer to vector data.
168: - For default PETSc vectors, VecGetArray() returns a pointer to
169: the data array. Otherwise, the routine is implementation dependent.
170: - You MUST call VecRestoreArray() when you no longer need access to
171: the array.
172: - Note that the Fortran interface to VecGetArray() differs from the
173: C version. See the users manual for details.
174: */
175: DMDAVecGetArray(da,U,&u);
177: /*
178: We initialize the solution array by simply writing the solution
179: directly into the array locations. Alternatively, we could use
180: VecSetValues() or VecSetValuesLocal().
181: */
182: for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
184: /*
185: Restore vector
186: */
187: DMDAVecRestoreArray(da,U,&u);
188: return 0;
189: }
190: /* --------------------------------------------------------------------- */
191: /*
192: Solution - Computes the exact solution at a given time.
194: Input Parameters:
195: t - current time
196: solution - vector in which exact solution will be computed
197: appctx - user-defined application context
199: Output Parameter:
200: solution - vector with the newly computed exact solution
201: */
202: PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx)
203: {
204: PetscScalar *u,ex1,ex2,sc1,sc2,h;
205: PetscInt i,mstart,mend,xm,M;
206: DM da;
208: TSGetDM(ts,&da);
209: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
210: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
211: h = 1.0/M;
212: mend = mstart + xm;
213: /*
214: Get a pointer to vector data.
215: */
216: DMDAVecGetArray(da,U,&u);
218: /*
219: Simply write the solution directly into the array locations.
220: Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
221: */
222: ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*appctx->d*t);
223: ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*appctx->d*t);
224: sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h;
225: 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;
227: /*
228: Restore vector
229: */
230: DMDAVecRestoreArray(da,U,&u);
231: return 0;
232: }
234: /* --------------------------------------------------------------------- */
235: /*
236: RHSMatrixHeat - User-provided routine to compute the right-hand-side
237: matrix for the heat equation.
239: Input Parameters:
240: ts - the TS context
241: t - current time
242: global_in - global input vector
243: dummy - optional user-defined context, as set by TSetRHSJacobian()
245: Output Parameters:
246: AA - Jacobian matrix
247: BB - optionally different preconditioning matrix
248: str - flag indicating matrix structure
250: Notes:
251: Recall that MatSetValues() uses 0-based row and column numbers
252: in Fortran as well as in C.
253: */
254: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec U,Mat AA,Mat BB,void *ctx)
255: {
256: Mat A = AA; /* Jacobian matrix */
257: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
258: PetscInt mstart, mend;
259: PetscInt i,idx[3],M,xm;
260: PetscScalar v[3],h;
261: DM da;
263: TSGetDM(ts,&da);
264: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
265: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
266: h = 1.0/M;
267: mend = mstart + xm;
268: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
269: Compute entries for the locally owned part of the matrix
270: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
271: /*
272: Set matrix rows corresponding to boundary data
273: */
275: /* diffusion */
276: v[0] = appctx->d/(h*h);
277: v[1] = -2.0*appctx->d/(h*h);
278: v[2] = appctx->d/(h*h);
279: if (!mstart) {
280: idx[0] = M-1; idx[1] = 0; idx[2] = 1;
281: MatSetValues(A,1,&mstart,3,idx,v,INSERT_VALUES);
282: mstart++;
283: }
285: if (mend == M) {
286: mend--;
287: idx[0] = M-2; idx[1] = M-1; idx[2] = 0;
288: MatSetValues(A,1,&mend,3,idx,v,INSERT_VALUES);
289: }
291: /*
292: Set matrix rows corresponding to interior data. We construct the
293: matrix one row at a time.
294: */
295: for (i=mstart; i<mend; i++) {
296: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
297: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
298: }
299: MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY);
300: MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY);
302: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
303: mend = mstart + xm;
304: if (!appctx->upwind) {
305: /* advection -- centered differencing */
306: v[0] = -.5*appctx->a/(h);
307: v[1] = .5*appctx->a/(h);
308: if (!mstart) {
309: idx[0] = M-1; idx[1] = 1;
310: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
311: mstart++;
312: }
314: if (mend == M) {
315: mend--;
316: idx[0] = M-2; idx[1] = 0;
317: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
318: }
320: for (i=mstart; i<mend; i++) {
321: idx[0] = i-1; idx[1] = i+1;
322: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
323: }
324: } else {
325: /* advection -- upwinding */
326: v[0] = -appctx->a/(h);
327: v[1] = appctx->a/(h);
328: if (!mstart) {
329: idx[0] = 0; idx[1] = 1;
330: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
331: mstart++;
332: }
334: if (mend == M) {
335: mend--;
336: idx[0] = M-1; idx[1] = 0;
337: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
338: }
340: for (i=mstart; i<mend; i++) {
341: idx[0] = i; idx[1] = i+1;
342: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
343: }
344: }
346: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
347: Complete the matrix assembly process and set some options
348: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
349: /*
350: Assemble matrix, using the 2-step process:
351: MatAssemblyBegin(), MatAssemblyEnd()
352: Computations can be done while messages are in transition
353: by placing code between these two statements.
354: */
355: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
356: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
358: /*
359: Set and option to indicate that we will never add a new nonzero location
360: to the matrix. If we do, it will generate an error.
361: */
362: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
363: return 0;
364: }
366: /*TEST
368: test:
369: args: -pc_type mg -da_refine 2 -ts_view -ts_monitor -ts_max_time .3 -mg_levels_ksp_max_it 3
370: requires: double
371: filter: grep -v "total number of"
373: test:
374: suffix: 2
375: 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
376: requires: x
377: output_file: output/ex3_1.out
378: requires: double
379: filter: grep -v "total number of"
381: TEST*/