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