Actual source code: ex3.c
petsc-3.12.2 2019-11-22
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);
128: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
129: Free work space. All PETSc objects should be destroyed when they
130: are no longer needed.
131: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
133: TSDestroy(&ts);
134: VecDestroy(&U);
135: DMDestroy(&da);
137: /*
138: Always call PetscFinalize() before exiting a program. This routine
139: - finalizes the PETSc libraries as well as MPI
140: - provides summary and diagnostic information if certain runtime
141: options are chosen (e.g., -log_view).
142: */
143: PetscFinalize();
144: return ierr;
145: }
146: /* --------------------------------------------------------------------- */
147: /*
148: InitialConditions - Computes the solution at the initial time.
150: Input Parameter:
151: u - uninitialized solution vector (global)
152: appctx - user-defined application context
154: Output Parameter:
155: u - vector with solution at initial time (global)
156: */
157: PetscErrorCode InitialConditions(TS ts,Vec U,AppCtx *appctx)
158: {
159: PetscScalar *u,h;
161: PetscInt i,mstart,mend,xm,M;
162: DM da;
164: TSGetDM(ts,&da);
165: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
166: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
167: h = 1.0/M;
168: mend = mstart + xm;
169: /*
170: Get a pointer to vector data.
171: - For default PETSc vectors, VecGetArray() returns a pointer to
172: the data array. Otherwise, the routine is implementation dependent.
173: - You MUST call VecRestoreArray() when you no longer need access to
174: the array.
175: - Note that the Fortran interface to VecGetArray() differs from the
176: C version. See the users manual for details.
177: */
178: DMDAVecGetArray(da,U,&u);
180: /*
181: We initialize the solution array by simply writing the solution
182: directly into the array locations. Alternatively, we could use
183: VecSetValues() or VecSetValuesLocal().
184: */
185: for (i=mstart; i<mend; i++) u[i] = PetscSinScalar(PETSC_PI*i*6.*h) + 3.*PetscSinScalar(PETSC_PI*i*2.*h);
187: /*
188: Restore vector
189: */
190: DMDAVecRestoreArray(da,U,&u);
191: return 0;
192: }
193: /* --------------------------------------------------------------------- */
194: /*
195: Solution - Computes the exact solution at a given time.
197: Input Parameters:
198: t - current time
199: solution - vector in which exact solution will be computed
200: appctx - user-defined application context
202: Output Parameter:
203: solution - vector with the newly computed exact solution
204: */
205: PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *appctx)
206: {
207: PetscScalar *u,ex1,ex2,sc1,sc2,h;
209: PetscInt i,mstart,mend,xm,M;
210: DM da;
212: TSGetDM(ts,&da);
213: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
214: DMDAGetInfo(da,PETSC_IGNORE,&M,0,0,0,0,0,0,0,0,0,0,0);
215: h = 1.0/M;
216: mend = mstart + xm;
217: /*
218: Get a pointer to vector data.
219: */
220: DMDAVecGetArray(da,U,&u);
222: /*
223: Simply write the solution directly into the array locations.
224: Alternatively, we culd use VecSetValues() or VecSetValuesLocal().
225: */
226: ex1 = PetscExpScalar(-36.*PETSC_PI*PETSC_PI*appctx->d*t);
227: ex2 = PetscExpScalar(-4.*PETSC_PI*PETSC_PI*appctx->d*t);
228: sc1 = PETSC_PI*6.*h; sc2 = PETSC_PI*2.*h;
229: 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;
231: /*
232: Restore vector
233: */
234: DMDAVecRestoreArray(da,U,&u);
235: return 0;
236: }
238: /* --------------------------------------------------------------------- */
239: /*
240: RHSMatrixHeat - User-provided routine to compute the right-hand-side
241: matrix for the heat equation.
243: Input Parameters:
244: ts - the TS context
245: t - current time
246: global_in - global input vector
247: dummy - optional user-defined context, as set by TSetRHSJacobian()
249: Output Parameters:
250: AA - Jacobian matrix
251: BB - optionally different preconditioning matrix
252: str - flag indicating matrix structure
254: Notes:
255: Recall that MatSetValues() uses 0-based row and column numbers
256: in Fortran as well as in C.
257: */
258: PetscErrorCode RHSMatrixHeat(TS ts,PetscReal t,Vec U,Mat AA,Mat BB,void *ctx)
259: {
260: Mat A = AA; /* Jacobian matrix */
261: AppCtx *appctx = (AppCtx*)ctx; /* user-defined application context */
262: PetscInt mstart, mend;
264: PetscInt i,idx[3],M,xm;
265: PetscScalar v[3],h;
266: DM da;
268: TSGetDM(ts,&da);
269: DMDAGetInfo(da,0,&M,0,0,0,0,0,0,0,0,0,0,0);
270: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
271: h = 1.0/M;
272: mend = mstart + xm;
273: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
274: Compute entries for the locally owned part of the matrix
275: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
276: /*
277: Set matrix rows corresponding to boundary data
278: */
280: /* diffusion */
281: v[0] = appctx->d/(h*h);
282: v[1] = -2.0*appctx->d/(h*h);
283: v[2] = appctx->d/(h*h);
284: if (!mstart) {
285: idx[0] = M-1; idx[1] = 0; idx[2] = 1;
286: MatSetValues(A,1,&mstart,3,idx,v,INSERT_VALUES);
287: mstart++;
288: }
290: if (mend == M) {
291: mend--;
292: idx[0] = M-2; idx[1] = M-1; idx[2] = 0;
293: MatSetValues(A,1,&mend,3,idx,v,INSERT_VALUES);
294: }
296: /*
297: Set matrix rows corresponding to interior data. We construct the
298: matrix one row at a time.
299: */
300: for (i=mstart; i<mend; i++) {
301: idx[0] = i-1; idx[1] = i; idx[2] = i+1;
302: MatSetValues(A,1,&i,3,idx,v,INSERT_VALUES);
303: }
304: MatAssemblyBegin(A,MAT_FLUSH_ASSEMBLY);
305: MatAssemblyEnd(A,MAT_FLUSH_ASSEMBLY);
307: DMDAGetCorners(da,&mstart,0,0,&xm,0,0);
308: mend = mstart + xm;
309: if (!appctx->upwind) {
310: /* advection -- centered differencing */
311: v[0] = -.5*appctx->a/(h);
312: v[1] = .5*appctx->a/(h);
313: if (!mstart) {
314: idx[0] = M-1; idx[1] = 1;
315: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
316: mstart++;
317: }
319: if (mend == M) {
320: mend--;
321: idx[0] = M-2; idx[1] = 0;
322: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
323: }
325: for (i=mstart; i<mend; i++) {
326: idx[0] = i-1; idx[1] = i+1;
327: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
328: }
329: } else {
330: /* advection -- upwinding */
331: v[0] = -appctx->a/(h);
332: v[1] = appctx->a/(h);
333: if (!mstart) {
334: idx[0] = 0; idx[1] = 1;
335: MatSetValues(A,1,&mstart,2,idx,v,ADD_VALUES);
336: mstart++;
337: }
339: if (mend == M) {
340: mend--;
341: idx[0] = M-1; idx[1] = 0;
342: MatSetValues(A,1,&mend,2,idx,v,ADD_VALUES);
343: }
345: for (i=mstart; i<mend; i++) {
346: idx[0] = i; idx[1] = i+1;
347: MatSetValues(A,1,&i,2,idx,v,ADD_VALUES);
348: }
349: }
352: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
353: Complete the matrix assembly process and set some options
354: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
355: /*
356: Assemble matrix, using the 2-step process:
357: MatAssemblyBegin(), MatAssemblyEnd()
358: Computations can be done while messages are in transition
359: by placing code between these two statements.
360: */
361: MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
362: MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
364: /*
365: Set and option to indicate that we will never add a new nonzero location
366: to the matrix. If we do, it will generate an error.
367: */
368: MatSetOption(A,MAT_NEW_NONZERO_LOCATION_ERR,PETSC_TRUE);
369: return 0;
370: }
373: /*TEST
375: test:
376: args: -pc_type mg -da_refine 2 -ts_view -ts_monitor -ts_max_time .3 -mg_levels_ksp_max_it 3
377: requires: double
379: test:
380: suffix: 2
381: 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
382: requires: x
383: output_file: output/ex3_1.out
384: requires: double
386: TEST*/