Actual source code: ex1.c
1: static char help[] = "Nonlinear Reaction Problem from Chemistry.\n";
3: /*F
5: This directory contains examples based on the PDES/ODES given in the book
6: Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
7: W. Hundsdorf and J.G. Verwer
9: Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry
11: \begin{eqnarray}
12: {U_1}_t - k U_1 U_2 & = & 0 \\
13: {U_2}_t - k U_1 U_2 & = & 0 \\
14: {U_3}_t - k U_1 U_2 & = & 0
15: \end{eqnarray}
17: Helpful runtime monitoring options:
18: -ts_view - prints information about the solver being used
19: -ts_monitor - prints the progress of the solver
20: -ts_adapt_monitor - prints the progress of the time-step adaptor
21: -ts_monitor_lg_timestep - plots the size of each timestep (at each time-step)
22: -ts_monitor_lg_solution - plots each component of the solution as a function of time (at each timestep)
23: -ts_monitor_lg_error - plots each component of the error in the solution as a function of time (at each timestep)
24: -draw_pause -2 - hold the plots a the end of the solution process, enter a mouse press in each window to end the process
26: -ts_monitor_lg_timestep -1 - plots the size of each timestep (at the end of the solution process)
27: -ts_monitor_lg_solution -1 - plots each component of the solution as a function of time (at the end of the solution process)
28: -ts_monitor_lg_error -1 - plots each component of the error in the solution as a function of time (at the end of the solution process)
29: -lg_use_markers false - do NOT show the data points on the plots
30: -draw_save - save the timestep and solution plot as a .Gif image file
32: F*/
34: /*
35: Project: Generate a nicely formatted HTML page using
36: 1) the HTML version of the source code and text in this file, use make html to generate the file ex1.c.html
37: 2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_$_1_0.Gif) and
38: 3) the text output (output.txt) generated by running the following commands.
39: 4) <iframe src="generated_topics.html" scrolling="no" frameborder="0" width=600 height=300></iframe>
41: rm -rf *.Gif
42: ./ex1 -ts_monitor_lg_error -1 -ts_monitor_lg_solution -1 -draw_pause -2 -ts_max_steps 100 -ts_monitor_lg_timestep -1 -draw_save -draw_virtual -ts_monitor -ts_adapt_monitor -ts_view > output.txt
44: For example something like
45: <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
46: <html>
47: <head>
48: <meta http-equiv="content-type" content="text/html;charset=utf-8">
49: <title>PETSc Example -- Nonlinear Reaction Problem from Chemistry</title>
50: </head>
51: <body>
52: <iframe src="ex1.c.html" scrolling="yes" frameborder="1" width=2000 height=400></iframe>
53: <img alt="" src="Draw_0x84000000_0_0.Gif"/><img alt="" src="Draw_0x84000001_0_0.Gif"/><img alt="" src="Draw_0x84000001_1_0.Gif"/>
54: <iframe src="output.txt" scrolling="yes" frameborder="1" width=2000 height=1000></iframe>
55: </body>
56: </html>
58: */
60: /*
61: Include "petscts.h" so that we can use TS solvers. Note that this
62: file automatically includes:
63: petscsys.h - base PETSc routines petscvec.h - vectors
64: petscmat.h - matrices
65: petscis.h - index sets petscksp.h - Krylov subspace methods
66: petscviewer.h - viewers petscpc.h - preconditioners
67: petscksp.h - linear solvers
68: */
70: #include <petscts.h>
72: typedef struct {
73: PetscScalar k;
74: Vec initialsolution;
75: } AppCtx;
77: PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v)
78: {
79: PetscFunctionBegin;
80: PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR));
81: PetscFunctionReturn(PETSC_SUCCESS);
82: }
84: PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v)
85: {
86: PetscFunctionBegin;
87: PetscCall(PetscNew(ctx));
88: PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR));
89: PetscFunctionReturn(PETSC_SUCCESS);
90: }
92: /*
93: Defines the ODE passed to the ODE solver
94: */
95: PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
96: {
97: PetscScalar *f;
98: const PetscScalar *u, *udot;
100: PetscFunctionBegin;
101: /* The next three lines allow us to access the entries of the vectors directly */
102: PetscCall(VecGetArrayRead(U, &u));
103: PetscCall(VecGetArrayRead(Udot, &udot));
104: PetscCall(VecGetArrayWrite(F, &f));
105: f[0] = udot[0] + ctx->k * u[0] * u[1];
106: f[1] = udot[1] + ctx->k * u[0] * u[1];
107: f[2] = udot[2] - ctx->k * u[0] * u[1];
108: PetscCall(VecRestoreArrayRead(U, &u));
109: PetscCall(VecRestoreArrayRead(Udot, &udot));
110: PetscCall(VecRestoreArrayWrite(F, &f));
111: PetscFunctionReturn(PETSC_SUCCESS);
112: }
114: /*
115: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
116: */
117: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
118: {
119: PetscInt rowcol[] = {0, 1, 2};
120: PetscScalar J[3][3];
121: const PetscScalar *u, *udot;
123: PetscFunctionBegin;
124: PetscCall(VecGetArrayRead(U, &u));
125: PetscCall(VecGetArrayRead(Udot, &udot));
126: J[0][0] = a + ctx->k * u[1];
127: J[0][1] = ctx->k * u[0];
128: J[0][2] = 0.0;
129: J[1][0] = ctx->k * u[1];
130: J[1][1] = a + ctx->k * u[0];
131: J[1][2] = 0.0;
132: J[2][0] = -ctx->k * u[1];
133: J[2][1] = -ctx->k * u[0];
134: J[2][2] = a;
135: PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
136: PetscCall(VecRestoreArrayRead(U, &u));
137: PetscCall(VecRestoreArrayRead(Udot, &udot));
139: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
140: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
141: if (A != B) {
142: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
143: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
144: }
145: PetscFunctionReturn(PETSC_SUCCESS);
146: }
148: /*
149: Defines the exact (analytic) solution to the ODE
150: */
151: static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx)
152: {
153: const PetscScalar *uinit;
154: PetscScalar *u, d0, q;
156: PetscFunctionBegin;
157: PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit));
158: PetscCall(VecGetArrayWrite(U, &u));
159: d0 = uinit[0] - uinit[1];
160: if (d0 == 0.0) q = ctx->k * t;
161: else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0;
162: u[0] = uinit[0] / (1.0 + uinit[1] * q);
163: u[1] = u[0] - d0;
164: u[2] = uinit[1] + uinit[2] - u[1];
165: PetscCall(VecRestoreArrayWrite(U, &u));
166: PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: int main(int argc, char **argv)
171: {
172: TS ts; /* ODE integrator */
173: Vec U; /* solution will be stored here */
174: Mat A; /* Jacobian matrix */
175: PetscMPIInt size;
176: PetscInt n = 3;
177: AppCtx ctx;
178: PetscScalar *u;
179: const char *const names[] = {"U1", "U2", "U3", NULL};
180: PetscBool mf = PETSC_FALSE;
182: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183: Initialize program
184: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185: PetscFunctionBeginUser;
186: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
187: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
188: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
190: PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &mf, NULL));
192: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193: Create necessary matrix and vectors
194: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
195: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
196: PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
197: PetscCall(MatSetFromOptions(A));
198: PetscCall(MatSetUp(A));
200: PetscCall(MatCreateVecs(A, &U, NULL));
202: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203: Set runtime options
204: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205: ctx.k = .9;
206: PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL));
207: PetscCall(VecDuplicate(U, &ctx.initialsolution));
208: PetscCall(VecGetArrayWrite(ctx.initialsolution, &u));
209: u[0] = 1;
210: u[1] = .7;
211: u[2] = 0;
212: PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u));
213: PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL));
215: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216: Create timestepping solver context
217: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
218: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
219: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
220: PetscCall(TSSetType(ts, TSROSW));
221: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &ctx));
222: if (!mf) PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx));
223: else PetscCall(TSSetIJacobian(ts, NULL, NULL, (TSIJacobianFn *)IJacobian, &ctx));
224: PetscCall(TSSetSolutionFunction(ts, (TSSolutionFn *)Solution, &ctx));
226: {
227: DM dm;
228: void *ptr;
229: PetscCall(TSGetDM(ts, &dm));
230: PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr));
231: PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr));
232: PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad));
233: PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad));
234: }
236: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237: Set initial conditions
238: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
239: PetscCall(Solution(ts, 0, U, &ctx));
240: PetscCall(TSSetSolution(ts, U));
242: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243: Set solver options
244: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
245: PetscCall(TSSetTimeStep(ts, .001));
246: PetscCall(TSSetMaxSteps(ts, 1000));
247: PetscCall(TSSetMaxTime(ts, 20.0));
248: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
249: PetscCall(TSSetFromOptions(ts));
250: PetscCall(TSMonitorLGSetVariableNames(ts, names));
252: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253: Solve nonlinear system
254: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255: PetscCall(TSSolve(ts, U));
257: PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD));
259: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260: Free work space. All PETSc objects should be destroyed when they are no longer needed.
261: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
262: PetscCall(VecDestroy(&ctx.initialsolution));
263: PetscCall(MatDestroy(&A));
264: PetscCall(VecDestroy(&U));
265: PetscCall(TSDestroy(&ts));
267: PetscCall(PetscFinalize());
268: return 0;
269: }
271: /*TEST
273: test:
274: args: -ts_view
275: requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
277: test:
278: suffix: 2
279: args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view
280: requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
281: output_file: output/ex1_1.out
283: test:
284: requires: !single
285: suffix: 3
286: args: -ts_view -snes_mf_operator
287: requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
289: TEST*/