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