1 static char help[] = "Nonlinear Reaction Problem from Chemistry.\n";
2
3 /*F
4
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
8
9 Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry
10
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}
16
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
25
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
31
32 F*/
33
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>
40
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
43
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>
57
58 */
59
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 */
69
70 #include <petscts.h>
71
72 typedef struct {
73 PetscScalar k;
74 Vec initialsolution;
75 } AppCtx;
76
IFunctionView(AppCtx * ctx,PetscViewer v)77 PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v)
78 {
79 PetscFunctionBegin;
80 PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR));
81 PetscFunctionReturn(PETSC_SUCCESS);
82 }
83
IFunctionLoad(AppCtx ** ctx,PetscViewer v)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 }
91
92 /*
93 Defines the ODE passed to the ODE solver
94 */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx * ctx)95 PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
96 {
97 PetscScalar *f;
98 const PetscScalar *u, *udot;
99
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 }
113
114 /*
115 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
116 */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx * ctx)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;
122
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));
138
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 }
147
148 /*
149 Defines the exact (analytic) solution to the ODE
150 */
Solution(TS ts,PetscReal t,Vec U,AppCtx * ctx)151 static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx)
152 {
153 const PetscScalar *uinit;
154 PetscScalar *u, d0, q;
155
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 }
169
main(int argc,char ** argv)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;
181
182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183 Initialize program
184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185 PetscFunctionBeginUser;
186 PetscCall(PetscInitialize(&argc, &argv, NULL, 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");
189
190 PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &mf, NULL));
191
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));
199
200 PetscCall(MatCreateVecs(A, &U, NULL));
201
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));
214
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));
225
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 }
235
236 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237 Set initial conditions
238 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
239 PetscCall(Solution(ts, 0, U, &ctx));
240 PetscCall(TSSetSolution(ts, U));
241
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));
251
252 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253 Solve nonlinear system
254 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
255 PetscCall(TSSolve(ts, U));
256
257 PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD));
258
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));
266
267 PetscCall(PetscFinalize());
268 return 0;
269 }
270
271 /*TEST
272
273 test:
274 args: -ts_view
275 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
276
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
282
283 test:
284 requires: !single
285 suffix: 3
286 args: -ts_view -snes_mf_operator
287 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
288
289 TEST*/
290