1c4762a1bSJed Brown static char help[] = "Nonlinear Reaction Problem from Chemistry.\n";
2c4762a1bSJed Brown
3c4762a1bSJed Brown /*F
4c4762a1bSJed Brown
5c4762a1bSJed Brown This directory contains examples based on the PDES/ODES given in the book
6c4762a1bSJed Brown Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
7c4762a1bSJed Brown W. Hundsdorf and J.G. Verwer
8c4762a1bSJed Brown
9c4762a1bSJed Brown Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry
10c4762a1bSJed Brown
11c4762a1bSJed Brown \begin{eqnarray}
12c4762a1bSJed Brown {U_1}_t - k U_1 U_2 & = & 0 \\
13c4762a1bSJed Brown {U_2}_t - k U_1 U_2 & = & 0 \\
14c4762a1bSJed Brown {U_3}_t - k U_1 U_2 & = & 0
15c4762a1bSJed Brown \end{eqnarray}
16c4762a1bSJed Brown
17c4762a1bSJed Brown Helpful runtime monitoring options:
18c4762a1bSJed Brown -ts_view - prints information about the solver being used
19da81f932SPierre Jolivet -ts_monitor - prints the progress of the solver
20c4762a1bSJed Brown -ts_adapt_monitor - prints the progress of the time-step adaptor
21c4762a1bSJed Brown -ts_monitor_lg_timestep - plots the size of each timestep (at each time-step)
22c4762a1bSJed Brown -ts_monitor_lg_solution - plots each component of the solution as a function of time (at each timestep)
23c4762a1bSJed Brown -ts_monitor_lg_error - plots each component of the error in the solution as a function of time (at each timestep)
24c4762a1bSJed Brown -draw_pause -2 - hold the plots a the end of the solution process, enter a mouse press in each window to end the process
25c4762a1bSJed Brown
26c4762a1bSJed Brown -ts_monitor_lg_timestep -1 - plots the size of each timestep (at the end of the solution process)
27c4762a1bSJed Brown -ts_monitor_lg_solution -1 - plots each component of the solution as a function of time (at the end of the solution process)
28c4762a1bSJed Brown -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)
29c4762a1bSJed Brown -lg_use_markers false - do NOT show the data points on the plots
30c4762a1bSJed Brown -draw_save - save the timestep and solution plot as a .Gif image file
31c4762a1bSJed Brown
32c4762a1bSJed Brown F*/
33c4762a1bSJed Brown
34c4762a1bSJed Brown /*
3535cb6cd3SPierre Jolivet Project: Generate a nicely formatted HTML page using
36c4762a1bSJed Brown 1) the HTML version of the source code and text in this file, use make html to generate the file ex1.c.html
371baa6e33SBarry Smith 2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_$_1_0.Gif) and
38c4762a1bSJed Brown 3) the text output (output.txt) generated by running the following commands.
39c4762a1bSJed Brown 4) <iframe src="generated_topics.html" scrolling="no" frameborder="0" width=600 height=300></iframe>
40c4762a1bSJed Brown
41c4762a1bSJed Brown rm -rf *.Gif
42c4762a1bSJed Brown ./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
43c4762a1bSJed Brown
44c4762a1bSJed Brown For example something like
45c4762a1bSJed Brown <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
46c4762a1bSJed Brown <html>
47c4762a1bSJed Brown <head>
48c4762a1bSJed Brown <meta http-equiv="content-type" content="text/html;charset=utf-8">
49c4762a1bSJed Brown <title>PETSc Example -- Nonlinear Reaction Problem from Chemistry</title>
50c4762a1bSJed Brown </head>
51c4762a1bSJed Brown <body>
52c4762a1bSJed Brown <iframe src="ex1.c.html" scrolling="yes" frameborder="1" width=2000 height=400></iframe>
53c4762a1bSJed Brown <img alt="" src="Draw_0x84000000_0_0.Gif"/><img alt="" src="Draw_0x84000001_0_0.Gif"/><img alt="" src="Draw_0x84000001_1_0.Gif"/>
54c4762a1bSJed Brown <iframe src="output.txt" scrolling="yes" frameborder="1" width=2000 height=1000></iframe>
55c4762a1bSJed Brown </body>
56c4762a1bSJed Brown </html>
57c4762a1bSJed Brown
58c4762a1bSJed Brown */
59c4762a1bSJed Brown
60c4762a1bSJed Brown /*
61c4762a1bSJed Brown Include "petscts.h" so that we can use TS solvers. Note that this
62c4762a1bSJed Brown file automatically includes:
63c4762a1bSJed Brown petscsys.h - base PETSc routines petscvec.h - vectors
64c4762a1bSJed Brown petscmat.h - matrices
65c4762a1bSJed Brown petscis.h - index sets petscksp.h - Krylov subspace methods
66c4762a1bSJed Brown petscviewer.h - viewers petscpc.h - preconditioners
67c4762a1bSJed Brown petscksp.h - linear solvers
68c4762a1bSJed Brown */
69c4762a1bSJed Brown
70c4762a1bSJed Brown #include <petscts.h>
71c4762a1bSJed Brown
72c4762a1bSJed Brown typedef struct {
73c4762a1bSJed Brown PetscScalar k;
74c4762a1bSJed Brown Vec initialsolution;
75c4762a1bSJed Brown } AppCtx;
76c4762a1bSJed Brown
IFunctionView(AppCtx * ctx,PetscViewer v)77d71ae5a4SJacob Faibussowitsch PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v)
78d71ae5a4SJacob Faibussowitsch {
79c4762a1bSJed Brown PetscFunctionBegin;
809566063dSJacob Faibussowitsch PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR));
813ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
82c4762a1bSJed Brown }
83c4762a1bSJed Brown
IFunctionLoad(AppCtx ** ctx,PetscViewer v)84d71ae5a4SJacob Faibussowitsch PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v)
85d71ae5a4SJacob Faibussowitsch {
86c4762a1bSJed Brown PetscFunctionBegin;
879566063dSJacob Faibussowitsch PetscCall(PetscNew(ctx));
889566063dSJacob Faibussowitsch PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR));
893ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
90c4762a1bSJed Brown }
91c4762a1bSJed Brown
92c4762a1bSJed Brown /*
93c4762a1bSJed Brown Defines the ODE passed to the ODE solver
94c4762a1bSJed Brown */
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx * ctx)95d71ae5a4SJacob Faibussowitsch PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx)
96d71ae5a4SJacob Faibussowitsch {
97c4762a1bSJed Brown PetscScalar *f;
98c4762a1bSJed Brown const PetscScalar *u, *udot;
99c4762a1bSJed Brown
100c4762a1bSJed Brown PetscFunctionBegin;
101c4762a1bSJed Brown /* The next three lines allow us to access the entries of the vectors directly */
1029566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(U, &u));
1039566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Udot, &udot));
1049566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(F, &f));
105c4762a1bSJed Brown f[0] = udot[0] + ctx->k * u[0] * u[1];
106c4762a1bSJed Brown f[1] = udot[1] + ctx->k * u[0] * u[1];
107c4762a1bSJed Brown f[2] = udot[2] - ctx->k * u[0] * u[1];
1089566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(U, &u));
1099566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Udot, &udot));
1109566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(F, &f));
1113ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
112c4762a1bSJed Brown }
113c4762a1bSJed Brown
114c4762a1bSJed Brown /*
115c4762a1bSJed Brown Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
116c4762a1bSJed Brown */
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx * ctx)117d71ae5a4SJacob Faibussowitsch PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
118d71ae5a4SJacob Faibussowitsch {
119c4762a1bSJed Brown PetscInt rowcol[] = {0, 1, 2};
120c4762a1bSJed Brown PetscScalar J[3][3];
121c4762a1bSJed Brown const PetscScalar *u, *udot;
122c4762a1bSJed Brown
123c4762a1bSJed Brown PetscFunctionBegin;
1249566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(U, &u));
1259566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Udot, &udot));
1269371c9d4SSatish Balay J[0][0] = a + ctx->k * u[1];
1279371c9d4SSatish Balay J[0][1] = ctx->k * u[0];
1289371c9d4SSatish Balay J[0][2] = 0.0;
1299371c9d4SSatish Balay J[1][0] = ctx->k * u[1];
1309371c9d4SSatish Balay J[1][1] = a + ctx->k * u[0];
1319371c9d4SSatish Balay J[1][2] = 0.0;
1329371c9d4SSatish Balay J[2][0] = -ctx->k * u[1];
1339371c9d4SSatish Balay J[2][1] = -ctx->k * u[0];
1349371c9d4SSatish Balay J[2][2] = a;
1359566063dSJacob Faibussowitsch PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
1369566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(U, &u));
1379566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Udot, &udot));
138c4762a1bSJed Brown
1399566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1409566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
141c4762a1bSJed Brown if (A != B) {
1429566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1439566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
144c4762a1bSJed Brown }
1453ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
146c4762a1bSJed Brown }
147c4762a1bSJed Brown
148c4762a1bSJed Brown /*
149c4762a1bSJed Brown Defines the exact (analytic) solution to the ODE
150c4762a1bSJed Brown */
Solution(TS ts,PetscReal t,Vec U,AppCtx * ctx)151d71ae5a4SJacob Faibussowitsch static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx)
152d71ae5a4SJacob Faibussowitsch {
153c4762a1bSJed Brown const PetscScalar *uinit;
154c4762a1bSJed Brown PetscScalar *u, d0, q;
155c4762a1bSJed Brown
156c4762a1bSJed Brown PetscFunctionBegin;
1579566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit));
1589566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(U, &u));
159c4762a1bSJed Brown d0 = uinit[0] - uinit[1];
160c4762a1bSJed Brown if (d0 == 0.0) q = ctx->k * t;
161c4762a1bSJed Brown else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0;
162c4762a1bSJed Brown u[0] = uinit[0] / (1.0 + uinit[1] * q);
163c4762a1bSJed Brown u[1] = u[0] - d0;
164c4762a1bSJed Brown u[2] = uinit[1] + uinit[2] - u[1];
1659566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(U, &u));
1669566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit));
1673ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS);
168c4762a1bSJed Brown }
169c4762a1bSJed Brown
main(int argc,char ** argv)170d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv)
171d71ae5a4SJacob Faibussowitsch {
172c4762a1bSJed Brown TS ts; /* ODE integrator */
173c4762a1bSJed Brown Vec U; /* solution will be stored here */
174c4762a1bSJed Brown Mat A; /* Jacobian matrix */
175c4762a1bSJed Brown PetscMPIInt size;
176c4762a1bSJed Brown PetscInt n = 3;
177c4762a1bSJed Brown AppCtx ctx;
178c4762a1bSJed Brown PetscScalar *u;
179c4762a1bSJed Brown const char *const names[] = {"U1", "U2", "U3", NULL};
18073a84a35SBarry Smith PetscBool mf = PETSC_FALSE;
181c4762a1bSJed Brown
182c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
183c4762a1bSJed Brown Initialize program
184c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
185327415f7SBarry Smith PetscFunctionBeginUser;
186*c8025a54SPierre Jolivet PetscCall(PetscInitialize(&argc, &argv, NULL, help));
1879566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1883c633725SBarry Smith PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
189c4762a1bSJed Brown
19073a84a35SBarry Smith PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &mf, NULL));
19173a84a35SBarry Smith
192c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
193c4762a1bSJed Brown Create necessary matrix and vectors
194c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1959566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
1969566063dSJacob Faibussowitsch PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
1979566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(A));
1989566063dSJacob Faibussowitsch PetscCall(MatSetUp(A));
199c4762a1bSJed Brown
2009566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A, &U, NULL));
201c4762a1bSJed Brown
202c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
203c4762a1bSJed Brown Set runtime options
204c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
205c4762a1bSJed Brown ctx.k = .9;
2069566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL));
2079566063dSJacob Faibussowitsch PetscCall(VecDuplicate(U, &ctx.initialsolution));
2089566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(ctx.initialsolution, &u));
209c4762a1bSJed Brown u[0] = 1;
210c4762a1bSJed Brown u[1] = .7;
211c4762a1bSJed Brown u[2] = 0;
2129566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u));
2139566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL));
214c4762a1bSJed Brown
215c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
216c4762a1bSJed Brown Create timestepping solver context
217c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2189566063dSJacob Faibussowitsch PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
2199566063dSJacob Faibussowitsch PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
2209566063dSJacob Faibussowitsch PetscCall(TSSetType(ts, TSROSW));
2218434afd1SBarry Smith PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &ctx));
22273a84a35SBarry Smith if (!mf) PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx));
22373a84a35SBarry Smith else PetscCall(TSSetIJacobian(ts, NULL, NULL, (TSIJacobianFn *)IJacobian, &ctx));
2248434afd1SBarry Smith PetscCall(TSSetSolutionFunction(ts, (TSSolutionFn *)Solution, &ctx));
225c4762a1bSJed Brown
226c4762a1bSJed Brown {
227c4762a1bSJed Brown DM dm;
228c4762a1bSJed Brown void *ptr;
2299566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &dm));
2309566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr));
2319566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr));
2329566063dSJacob Faibussowitsch PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode (*)(void *, PetscViewer))IFunctionView, (PetscErrorCode (*)(void **, PetscViewer))IFunctionLoad));
2339566063dSJacob Faibussowitsch PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode (*)(void *, PetscViewer))IFunctionView, (PetscErrorCode (*)(void **, PetscViewer))IFunctionLoad));
234c4762a1bSJed Brown }
235c4762a1bSJed Brown
236c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
237c4762a1bSJed Brown Set initial conditions
238c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2399566063dSJacob Faibussowitsch PetscCall(Solution(ts, 0, U, &ctx));
2409566063dSJacob Faibussowitsch PetscCall(TSSetSolution(ts, U));
241c4762a1bSJed Brown
242c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
243c4762a1bSJed Brown Set solver options
244c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2459566063dSJacob Faibussowitsch PetscCall(TSSetTimeStep(ts, .001));
2469566063dSJacob Faibussowitsch PetscCall(TSSetMaxSteps(ts, 1000));
2479566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts, 20.0));
2489566063dSJacob Faibussowitsch PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
2499566063dSJacob Faibussowitsch PetscCall(TSSetFromOptions(ts));
2509566063dSJacob Faibussowitsch PetscCall(TSMonitorLGSetVariableNames(ts, names));
251c4762a1bSJed Brown
252c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
253c4762a1bSJed Brown Solve nonlinear system
254c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2559566063dSJacob Faibussowitsch PetscCall(TSSolve(ts, U));
256c4762a1bSJed Brown
2579566063dSJacob Faibussowitsch PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD));
258c4762a1bSJed Brown
259c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
260c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they are no longer needed.
261c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2629566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ctx.initialsolution));
2639566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A));
2649566063dSJacob Faibussowitsch PetscCall(VecDestroy(&U));
2659566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts));
266c4762a1bSJed Brown
2679566063dSJacob Faibussowitsch PetscCall(PetscFinalize());
268b122ec5aSJacob Faibussowitsch return 0;
269c4762a1bSJed Brown }
270c4762a1bSJed Brown
271c4762a1bSJed Brown /*TEST
272c4762a1bSJed Brown
273c4762a1bSJed Brown test:
274c4762a1bSJed Brown args: -ts_view
275dfd57a17SPierre Jolivet requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
276c4762a1bSJed Brown
277c4762a1bSJed Brown test:
278c4762a1bSJed Brown suffix: 2
279c4762a1bSJed Brown args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view
280dfd57a17SPierre Jolivet requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
281c4762a1bSJed Brown output_file: output/ex1_1.out
282c4762a1bSJed Brown
28373a84a35SBarry Smith test:
28473a84a35SBarry Smith requires: !single
28573a84a35SBarry Smith suffix: 3
28673a84a35SBarry Smith args: -ts_view -snes_mf_operator
28773a84a35SBarry Smith requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
28873a84a35SBarry Smith
289c4762a1bSJed Brown TEST*/
290