Actual source code: ex12.c
2: static char help[] = "Tests PetscObjectSetOptions() for TS object\n\n";
4: /* ------------------------------------------------------------------------
6: This program solves the PDE
8: u * u_xx
9: u_t = ---------
10: 2*(t+1)^2
12: on the domain 0 <= x <= 1, with boundary conditions
13: u(t,0) = t + 1, u(t,1) = 2*t + 2,
14: and initial condition
15: u(0,x) = 1 + x*x.
17: The exact solution is:
18: u(t,x) = (1 + x*x) * (1 + t)
20: Note that since the solution is linear in time and quadratic in x,
21: the finite difference scheme actually computes the "exact" solution.
23: We use by default the backward Euler method.
25: ------------------------------------------------------------------------- */
27: /*
28: Include "petscts.h" to use the PETSc timestepping routines. Note that
29: this file automatically includes "petscsys.h" and other lower-level
30: PETSc include files.
32: Include the "petscdmda.h" to allow us to use the distributed array data
33: structures to manage the parallel grid.
34: */
35: #include <petscts.h>
36: #include <petscdm.h>
37: #include <petscdmda.h>
38: #include <petscdraw.h>
40: /*
41: User-defined application context - contains data needed by the
42: application-provided callback routines.
43: */
44: typedef struct {
45: MPI_Comm comm; /* communicator */
46: DM da; /* distributed array data structure */
47: Vec localwork; /* local ghosted work vector */
48: Vec u_local; /* local ghosted approximate solution vector */
49: Vec solution; /* global exact solution vector */
50: PetscInt m; /* total number of grid points */
51: PetscReal h; /* mesh width: h = 1/(m-1) */
52: } AppCtx;
54: /*
55: User-defined routines, provided below.
56: */
57: extern PetscErrorCode InitialConditions(Vec, AppCtx *);
58: extern PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
59: extern PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
60: extern PetscErrorCode ExactSolution(PetscReal, Vec, AppCtx *);
62: int main(int argc, char **argv)
63: {
64: AppCtx appctx; /* user-defined application context */
65: TS ts; /* timestepping context */
66: Mat A; /* Jacobian matrix data structure */
67: Vec u; /* approximate solution vector */
68: PetscInt time_steps_max = 100; /* default max timesteps */
69: PetscReal dt;
70: PetscReal time_total_max = 100.0; /* default max total time */
71: PetscOptions options, optionscopy;
73: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
74: Initialize program and set problem parameters
75: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
77: PetscFunctionBeginUser;
78: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
80: PetscCall(PetscOptionsCreate(&options));
81: PetscCall(PetscOptionsSetValue(options, "-ts_monitor", "ascii"));
82: PetscCall(PetscOptionsSetValue(options, "-snes_monitor", "ascii"));
83: PetscCall(PetscOptionsSetValue(options, "-ksp_monitor", "ascii"));
85: appctx.comm = PETSC_COMM_WORLD;
86: appctx.m = 60;
88: PetscCall(PetscOptionsGetInt(options, NULL, "-M", &appctx.m, NULL));
90: appctx.h = 1.0 / (appctx.m - 1.0);
92: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
93: Create vector data structures
94: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
96: /*
97: Create distributed array (DMDA) to manage parallel grid and vectors
98: and to set up the ghost point communication pattern. There are M
99: total grid values spread equally among all the processors.
100: */
101: PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, appctx.m, 1, 1, NULL, &appctx.da));
102: PetscCall(PetscObjectSetOptions((PetscObject)appctx.da, options));
103: PetscCall(DMSetFromOptions(appctx.da));
104: PetscCall(DMSetUp(appctx.da));
106: /*
107: Extract global and local vectors from DMDA; we use these to store the
108: approximate solution. Then duplicate these for remaining vectors that
109: have the same types.
110: */
111: PetscCall(DMCreateGlobalVector(appctx.da, &u));
112: PetscCall(DMCreateLocalVector(appctx.da, &appctx.u_local));
114: /*
115: Create local work vector for use in evaluating right-hand-side function;
116: create global work vector for storing exact solution.
117: */
118: PetscCall(VecDuplicate(appctx.u_local, &appctx.localwork));
119: PetscCall(VecDuplicate(u, &appctx.solution));
121: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
122: Create timestepping solver context; set callback routine for
123: right-hand-side function evaluation.
124: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
127: PetscCall(PetscObjectSetOptions((PetscObject)ts, options));
128: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
129: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, &appctx));
131: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
132: For nonlinear problems, the user can provide a Jacobian evaluation
133: routine (or use a finite differencing approximation).
135: Create matrix data structure; set Jacobian evaluation routine.
136: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
138: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
139: PetscCall(PetscObjectSetOptions((PetscObject)A, options));
140: PetscCall(MatSetSizes(A, PETSC_DECIDE, PETSC_DECIDE, appctx.m, appctx.m));
141: PetscCall(MatSetFromOptions(A));
142: PetscCall(MatSetUp(A));
143: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, &appctx));
145: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
146: Set solution vector and initial timestep
147: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: dt = appctx.h / 2.0;
150: PetscCall(TSSetTimeStep(ts, dt));
152: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
153: Customize timestepping solver:
154: - Set the solution method to be the Backward Euler method.
155: - Set timestepping duration info
156: Then set runtime options, which can override these defaults.
157: For example,
158: -ts_max_steps <maxsteps> -ts_max_time <maxtime>
159: to override the defaults set by TSSetMaxSteps()/TSSetMaxTime().
160: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
162: PetscCall(TSSetType(ts, TSBEULER));
163: PetscCall(TSSetMaxSteps(ts, time_steps_max));
164: PetscCall(TSSetMaxTime(ts, time_total_max));
165: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
166: PetscCall(TSSetFromOptions(ts));
168: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
169: Solve the problem
170: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
172: /*
173: Evaluate initial conditions
174: */
175: PetscCall(InitialConditions(u, &appctx));
177: /*
178: Run the timestepping solver
179: */
180: PetscCall(TSSolve(ts, u));
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Free work space. All PETSc objects should be destroyed when they
184: are no longer needed.
185: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187: PetscCall(PetscObjectGetOptions((PetscObject)ts, &optionscopy));
188: PetscCheck(options == optionscopy, PETSC_COMM_WORLD, PETSC_ERR_PLIB, "PetscObjectGetOptions() failed");
190: PetscCall(TSDestroy(&ts));
191: PetscCall(VecDestroy(&u));
192: PetscCall(MatDestroy(&A));
193: PetscCall(DMDestroy(&appctx.da));
194: PetscCall(VecDestroy(&appctx.localwork));
195: PetscCall(VecDestroy(&appctx.solution));
196: PetscCall(VecDestroy(&appctx.u_local));
197: PetscCall(PetscOptionsDestroy(&options));
199: /*
200: Always call PetscFinalize() before exiting a program. This routine
201: - finalizes the PETSc libraries as well as MPI
202: - provides summary and diagnostic information if certain runtime
203: options are chosen (e.g., -log_view).
204: */
205: PetscCall(PetscFinalize());
206: return 0;
207: }
208: /* --------------------------------------------------------------------- */
209: /*
210: InitialConditions - Computes the solution at the initial time.
212: Input Parameters:
213: u - uninitialized solution vector (global)
214: appctx - user-defined application context
216: Output Parameter:
217: u - vector with solution at initial time (global)
218: */
219: PetscErrorCode InitialConditions(Vec u, AppCtx *appctx)
220: {
221: PetscScalar *u_localptr, h = appctx->h, x;
222: PetscInt i, mybase, myend;
224: PetscFunctionBeginUser;
225: /*
226: Determine starting point of each processor's range of
227: grid values.
228: */
229: PetscCall(VecGetOwnershipRange(u, &mybase, &myend));
231: /*
232: Get a pointer to vector data.
233: - For default PETSc vectors, VecGetArray() returns a pointer to
234: the data array. Otherwise, the routine is implementation dependent.
235: - You MUST call VecRestoreArray() when you no longer need access to
236: the array.
237: - Note that the Fortran interface to VecGetArray() differs from the
238: C version. See the users manual for details.
239: */
240: PetscCall(VecGetArray(u, &u_localptr));
242: /*
243: We initialize the solution array by simply writing the solution
244: directly into the array locations. Alternatively, we could use
245: VecSetValues() or VecSetValuesLocal().
246: */
247: for (i = mybase; i < myend; i++) {
248: x = h * (PetscReal)i; /* current location in global grid */
249: u_localptr[i - mybase] = 1.0 + x * x;
250: }
252: /*
253: Restore vector
254: */
255: PetscCall(VecRestoreArray(u, &u_localptr));
257: PetscFunctionReturn(PETSC_SUCCESS);
258: }
259: /* --------------------------------------------------------------------- */
260: /*
261: ExactSolution - Computes the exact solution at a given time.
263: Input Parameters:
264: t - current time
265: solution - vector in which exact solution will be computed
266: appctx - user-defined application context
268: Output Parameter:
269: solution - vector with the newly computed exact solution
270: */
271: PetscErrorCode ExactSolution(PetscReal t, Vec solution, AppCtx *appctx)
272: {
273: PetscScalar *s_localptr, h = appctx->h, x;
274: PetscInt i, mybase, myend;
276: PetscFunctionBeginUser;
277: /*
278: Determine starting and ending points of each processor's
279: range of grid values
280: */
281: PetscCall(VecGetOwnershipRange(solution, &mybase, &myend));
283: /*
284: Get a pointer to vector data.
285: */
286: PetscCall(VecGetArray(solution, &s_localptr));
288: /*
289: Simply write the solution directly into the array locations.
290: Alternatively, we could use VecSetValues() or VecSetValuesLocal().
291: */
292: for (i = mybase; i < myend; i++) {
293: x = h * (PetscReal)i;
294: s_localptr[i - mybase] = (t + 1.0) * (1.0 + x * x);
295: }
297: /*
298: Restore vector
299: */
300: PetscCall(VecRestoreArray(solution, &s_localptr));
301: PetscFunctionReturn(PETSC_SUCCESS);
302: }
303: /* --------------------------------------------------------------------- */
304: /*
305: RHSFunction - User-provided routine that evalues the right-hand-side
306: function of the ODE. This routine is set in the main program by
307: calling TSSetRHSFunction(). We compute:
308: global_out = F(global_in)
310: Input Parameters:
311: ts - timesteping context
312: t - current time
313: global_in - vector containing the current iterate
314: ctx - (optional) user-provided context for function evaluation.
315: In this case we use the appctx defined above.
317: Output Parameter:
318: global_out - vector containing the newly evaluated function
319: */
320: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec global_in, Vec global_out, void *ctx)
321: {
322: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
323: DM da = appctx->da; /* distributed array */
324: Vec local_in = appctx->u_local; /* local ghosted input vector */
325: Vec localwork = appctx->localwork; /* local ghosted work vector */
326: PetscInt i, localsize;
327: PetscMPIInt rank, size;
328: PetscScalar *copyptr, sc;
329: const PetscScalar *localptr;
331: PetscFunctionBeginUser;
332: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
333: Get ready for local function computations
334: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
335: /*
336: Scatter ghost points to local vector, using the 2-step process
337: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
338: By placing code between these two statements, computations can be
339: done while messages are in transition.
340: */
341: PetscCall(DMGlobalToLocalBegin(da, global_in, INSERT_VALUES, local_in));
342: PetscCall(DMGlobalToLocalEnd(da, global_in, INSERT_VALUES, local_in));
344: /*
345: Access directly the values in our local INPUT work array
346: */
347: PetscCall(VecGetArrayRead(local_in, &localptr));
349: /*
350: Access directly the values in our local OUTPUT work array
351: */
352: PetscCall(VecGetArray(localwork, ©ptr));
354: sc = 1.0 / (appctx->h * appctx->h * 2.0 * (1.0 + t) * (1.0 + t));
356: /*
357: Evaluate our function on the nodes owned by this processor
358: */
359: PetscCall(VecGetLocalSize(local_in, &localsize));
361: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
362: Compute entries for the locally owned part
363: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
365: /*
366: Handle boundary conditions: This is done by using the boundary condition
367: u(t,boundary) = g(t,boundary)
368: for some function g. Now take the derivative with respect to t to obtain
369: u_{t}(t,boundary) = g_{t}(t,boundary)
371: In our case, u(t,0) = t + 1, so that u_{t}(t,0) = 1
372: and u(t,1) = 2t+ 2, so that u_{t}(t,1) = 2
373: */
374: PetscCallMPI(MPI_Comm_rank(appctx->comm, &rank));
375: PetscCallMPI(MPI_Comm_size(appctx->comm, &size));
376: if (rank == 0) copyptr[0] = 1.0;
377: if (rank == size - 1) copyptr[localsize - 1] = 2.0;
379: /*
380: Handle the interior nodes where the PDE is replace by finite
381: difference operators.
382: */
383: for (i = 1; i < localsize - 1; i++) copyptr[i] = localptr[i] * sc * (localptr[i + 1] + localptr[i - 1] - 2.0 * localptr[i]);
385: /*
386: Restore vectors
387: */
388: PetscCall(VecRestoreArrayRead(local_in, &localptr));
389: PetscCall(VecRestoreArray(localwork, ©ptr));
391: /*
392: Insert values from the local OUTPUT vector into the global
393: output vector
394: */
395: PetscCall(DMLocalToGlobalBegin(da, localwork, INSERT_VALUES, global_out));
396: PetscCall(DMLocalToGlobalEnd(da, localwork, INSERT_VALUES, global_out));
398: PetscFunctionReturn(PETSC_SUCCESS);
399: }
400: /* --------------------------------------------------------------------- */
401: /*
402: RHSJacobian - User-provided routine to compute the Jacobian of
403: the nonlinear right-hand-side function of the ODE.
405: Input Parameters:
406: ts - the TS context
407: t - current time
408: global_in - global input vector
409: dummy - optional user-defined context, as set by TSetRHSJacobian()
411: Output Parameters:
412: AA - Jacobian matrix
413: BB - optionally different preconditioning matrix
414: str - flag indicating matrix structure
416: Notes:
417: RHSJacobian computes entries for the locally owned part of the Jacobian.
418: - Currently, all PETSc parallel matrix formats are partitioned by
419: contiguous chunks of rows across the processors.
420: - Each processor needs to insert only elements that it owns
421: locally (but any non-local elements will be sent to the
422: appropriate processor during matrix assembly).
423: - Always specify global row and columns of matrix entries when
424: using MatSetValues().
425: - Here, we set all entries for a particular row at once.
426: - Note that MatSetValues() uses 0-based row and column numbers
427: in Fortran as well as in C.
428: */
429: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec global_in, Mat AA, Mat BB, void *ctx)
430: {
431: AppCtx *appctx = (AppCtx *)ctx; /* user-defined application context */
432: Vec local_in = appctx->u_local; /* local ghosted input vector */
433: DM da = appctx->da; /* distributed array */
434: PetscScalar v[3], sc;
435: const PetscScalar *localptr;
436: PetscInt i, mstart, mend, mstarts, mends, idx[3], is;
438: PetscFunctionBeginUser;
439: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
440: Get ready for local Jacobian computations
441: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
442: /*
443: Scatter ghost points to local vector, using the 2-step process
444: DMGlobalToLocalBegin(), DMGlobalToLocalEnd().
445: By placing code between these two statements, computations can be
446: done while messages are in transition.
447: */
448: PetscCall(DMGlobalToLocalBegin(da, global_in, INSERT_VALUES, local_in));
449: PetscCall(DMGlobalToLocalEnd(da, global_in, INSERT_VALUES, local_in));
451: /*
452: Get pointer to vector data
453: */
454: PetscCall(VecGetArrayRead(local_in, &localptr));
456: /*
457: Get starting and ending locally owned rows of the matrix
458: */
459: PetscCall(MatGetOwnershipRange(BB, &mstarts, &mends));
460: mstart = mstarts;
461: mend = mends;
463: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
464: Compute entries for the locally owned part of the Jacobian.
465: - Currently, all PETSc parallel matrix formats are partitioned by
466: contiguous chunks of rows across the processors.
467: - Each processor needs to insert only elements that it owns
468: locally (but any non-local elements will be sent to the
469: appropriate processor during matrix assembly).
470: - Here, we set all entries for a particular row at once.
471: - We can set matrix entries either using either
472: MatSetValuesLocal() or MatSetValues().
473: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
475: /*
476: Set matrix rows corresponding to boundary data
477: */
478: if (mstart == 0) {
479: v[0] = 0.0;
480: PetscCall(MatSetValues(BB, 1, &mstart, 1, &mstart, v, INSERT_VALUES));
481: mstart++;
482: }
483: if (mend == appctx->m) {
484: mend--;
485: v[0] = 0.0;
486: PetscCall(MatSetValues(BB, 1, &mend, 1, &mend, v, INSERT_VALUES));
487: }
489: /*
490: Set matrix rows corresponding to interior data. We construct the
491: matrix one row at a time.
492: */
493: sc = 1.0 / (appctx->h * appctx->h * 2.0 * (1.0 + t) * (1.0 + t));
494: for (i = mstart; i < mend; i++) {
495: idx[0] = i - 1;
496: idx[1] = i;
497: idx[2] = i + 1;
498: is = i - mstart + 1;
499: v[0] = sc * localptr[is];
500: v[1] = sc * (localptr[is + 1] + localptr[is - 1] - 4.0 * localptr[is]);
501: v[2] = sc * localptr[is];
502: PetscCall(MatSetValues(BB, 1, &i, 3, idx, v, INSERT_VALUES));
503: }
505: /*
506: Restore vector
507: */
508: PetscCall(VecRestoreArrayRead(local_in, &localptr));
510: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
511: Complete the matrix assembly process and set some options
512: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
513: /*
514: Assemble matrix, using the 2-step process:
515: MatAssemblyBegin(), MatAssemblyEnd()
516: Computations can be done while messages are in transition
517: by placing code between these two statements.
518: */
519: PetscCall(MatAssemblyBegin(BB, MAT_FINAL_ASSEMBLY));
520: PetscCall(MatAssemblyEnd(BB, MAT_FINAL_ASSEMBLY));
521: if (BB != AA) {
522: PetscCall(MatAssemblyBegin(AA, MAT_FINAL_ASSEMBLY));
523: PetscCall(MatAssemblyEnd(AA, MAT_FINAL_ASSEMBLY));
524: }
526: /*
527: Set and option to indicate that we will never add a new nonzero location
528: to the matrix. If we do, it will generate an error.
529: */
530: PetscCall(MatSetOption(BB, MAT_NEW_NONZERO_LOCATION_ERR, PETSC_TRUE));
532: PetscFunctionReturn(PETSC_SUCCESS);
533: }
535: /*TEST
537: test:
538: requires: !single
540: TEST*/