xref: /petsc/src/ts/tutorials/advection-diffusion-reaction/ex1.c (revision 9371c9d470a9602b6d10a8bf50c9b2280a79e45a) 
1c4762a1bSJed Brown 
2c4762a1bSJed Brown static char help[] = "Nonlinear Reaction Problem from Chemistry.\n";
3c4762a1bSJed Brown 
4c4762a1bSJed Brown /*F
5c4762a1bSJed Brown 
6c4762a1bSJed Brown      This directory contains examples based on the PDES/ODES given in the book
7c4762a1bSJed Brown       Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
8c4762a1bSJed Brown       W. Hundsdorf and J.G. Verwer
9c4762a1bSJed Brown 
10c4762a1bSJed Brown      Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry
11c4762a1bSJed Brown 
12c4762a1bSJed Brown \begin{eqnarray}
13c4762a1bSJed Brown                  {U_1}_t  - k U_1 U_2  & = & 0 \\
14c4762a1bSJed Brown                  {U_2}_t  - k U_1 U_2 & = & 0 \\
15c4762a1bSJed Brown                  {U_3}_t  - k U_1 U_2 & = & 0
16c4762a1bSJed Brown \end{eqnarray}
17c4762a1bSJed Brown 
18c4762a1bSJed Brown      Helpful runtime monitoring options:
19c4762a1bSJed Brown          -ts_view                  -  prints information about the solver being used
20c4762a1bSJed Brown          -ts_monitor               -  prints the progess of the solver
21c4762a1bSJed Brown          -ts_adapt_monitor         -  prints the progress of the time-step adaptor
22c4762a1bSJed Brown          -ts_monitor_lg_timestep   -  plots the size of each timestep (at each time-step)
23c4762a1bSJed Brown          -ts_monitor_lg_solution   -  plots each component of the solution as a function of time (at each timestep)
24c4762a1bSJed Brown          -ts_monitor_lg_error      -  plots each component of the error in the solution as a function of time (at each timestep)
25c4762a1bSJed 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
26c4762a1bSJed Brown 
27c4762a1bSJed Brown          -ts_monitor_lg_timestep -1  -  plots the size of each timestep (at the end of the solution process)
28c4762a1bSJed Brown          -ts_monitor_lg_solution -1  -  plots each component of the solution as a function of time (at the end of the solution process)
29c4762a1bSJed 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)
30c4762a1bSJed Brown          -lg_use_markers false       -  do NOT show the data points on the plots
31c4762a1bSJed Brown          -draw_save                  -  save the timestep and solution plot as a .Gif image file
32c4762a1bSJed Brown 
33c4762a1bSJed Brown F*/
34c4762a1bSJed Brown 
35c4762a1bSJed Brown /*
36c4762a1bSJed Brown       Project: Generate a nicely formated HTML page using
37c4762a1bSJed Brown          1) the HTML version of the source code and text in this file, use make html to generate the file ex1.c.html
381baa6e33SBarry Smith          2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_$_1_0.Gif) and
39c4762a1bSJed Brown          3) the text output (output.txt) generated by running the following commands.
40c4762a1bSJed Brown          4) <iframe src="generated_topics.html" scrolling="no" frameborder="0"  width=600 height=300></iframe>
41c4762a1bSJed Brown 
42c4762a1bSJed Brown       rm -rf *.Gif
43c4762a1bSJed 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
44c4762a1bSJed Brown 
45c4762a1bSJed Brown       For example something like
46c4762a1bSJed Brown <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
47c4762a1bSJed Brown <html>
48c4762a1bSJed Brown   <head>
49c4762a1bSJed Brown     <meta http-equiv="content-type" content="text/html;charset=utf-8">
50c4762a1bSJed Brown     <title>PETSc Example -- Nonlinear Reaction Problem from Chemistry</title>
51c4762a1bSJed Brown   </head>
52c4762a1bSJed Brown   <body>
53c4762a1bSJed Brown   <iframe src="ex1.c.html" scrolling="yes" frameborder="1"  width=2000 height=400></iframe>
54c4762a1bSJed 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"/>
55c4762a1bSJed Brown   <iframe src="output.txt" scrolling="yes" frameborder="1"  width=2000 height=1000></iframe>
56c4762a1bSJed Brown   </body>
57c4762a1bSJed Brown </html>
58c4762a1bSJed Brown 
59c4762a1bSJed Brown */
60c4762a1bSJed Brown 
61c4762a1bSJed Brown /*
62c4762a1bSJed Brown    Include "petscts.h" so that we can use TS solvers.  Note that this
63c4762a1bSJed Brown    file automatically includes:
64c4762a1bSJed Brown      petscsys.h       - base PETSc routines   petscvec.h - vectors
65c4762a1bSJed Brown      petscmat.h - matrices
66c4762a1bSJed Brown      petscis.h     - index sets            petscksp.h - Krylov subspace methods
67c4762a1bSJed Brown      petscviewer.h - viewers               petscpc.h  - preconditioners
68c4762a1bSJed Brown      petscksp.h   - linear solvers
69c4762a1bSJed Brown */
70c4762a1bSJed Brown 
71c4762a1bSJed Brown #include <petscts.h>
72c4762a1bSJed Brown 
73c4762a1bSJed Brown typedef struct {
74c4762a1bSJed Brown   PetscScalar k;
75c4762a1bSJed Brown   Vec         initialsolution;
76c4762a1bSJed Brown } AppCtx;
77c4762a1bSJed Brown 
78*9371c9d4SSatish Balay PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v) {
79c4762a1bSJed Brown   PetscFunctionBegin;
809566063dSJacob Faibussowitsch   PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR));
81c4762a1bSJed Brown   PetscFunctionReturn(0);
82c4762a1bSJed Brown }
83c4762a1bSJed Brown 
84*9371c9d4SSatish Balay PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v) {
85c4762a1bSJed Brown   PetscFunctionBegin;
869566063dSJacob Faibussowitsch   PetscCall(PetscNew(ctx));
879566063dSJacob Faibussowitsch   PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR));
88c4762a1bSJed Brown   PetscFunctionReturn(0);
89c4762a1bSJed Brown }
90c4762a1bSJed Brown 
91c4762a1bSJed Brown /*
92c4762a1bSJed Brown      Defines the ODE passed to the ODE solver
93c4762a1bSJed Brown */
94*9371c9d4SSatish Balay PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx) {
95c4762a1bSJed Brown   PetscScalar       *f;
96c4762a1bSJed Brown   const PetscScalar *u, *udot;
97c4762a1bSJed Brown 
98c4762a1bSJed Brown   PetscFunctionBegin;
99c4762a1bSJed Brown   /*  The next three lines allow us to access the entries of the vectors directly */
1009566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(U, &u));
1019566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(Udot, &udot));
1029566063dSJacob Faibussowitsch   PetscCall(VecGetArrayWrite(F, &f));
103c4762a1bSJed Brown   f[0] = udot[0] + ctx->k * u[0] * u[1];
104c4762a1bSJed Brown   f[1] = udot[1] + ctx->k * u[0] * u[1];
105c4762a1bSJed Brown   f[2] = udot[2] - ctx->k * u[0] * u[1];
1069566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(U, &u));
1079566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(Udot, &udot));
1089566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayWrite(F, &f));
109c4762a1bSJed Brown   PetscFunctionReturn(0);
110c4762a1bSJed Brown }
111c4762a1bSJed Brown 
112c4762a1bSJed Brown /*
113c4762a1bSJed Brown      Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
114c4762a1bSJed Brown */
115*9371c9d4SSatish Balay PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx) {
116c4762a1bSJed Brown   PetscInt           rowcol[] = {0, 1, 2};
117c4762a1bSJed Brown   PetscScalar        J[3][3];
118c4762a1bSJed Brown   const PetscScalar *u, *udot;
119c4762a1bSJed Brown 
120c4762a1bSJed Brown   PetscFunctionBegin;
1219566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(U, &u));
1229566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(Udot, &udot));
123*9371c9d4SSatish Balay   J[0][0] = a + ctx->k * u[1];
124*9371c9d4SSatish Balay   J[0][1] = ctx->k * u[0];
125*9371c9d4SSatish Balay   J[0][2] = 0.0;
126*9371c9d4SSatish Balay   J[1][0] = ctx->k * u[1];
127*9371c9d4SSatish Balay   J[1][1] = a + ctx->k * u[0];
128*9371c9d4SSatish Balay   J[1][2] = 0.0;
129*9371c9d4SSatish Balay   J[2][0] = -ctx->k * u[1];
130*9371c9d4SSatish Balay   J[2][1] = -ctx->k * u[0];
131*9371c9d4SSatish Balay   J[2][2] = a;
1329566063dSJacob Faibussowitsch   PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES));
1339566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(U, &u));
1349566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(Udot, &udot));
135c4762a1bSJed Brown 
1369566063dSJacob Faibussowitsch   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
1379566063dSJacob Faibussowitsch   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
138c4762a1bSJed Brown   if (A != B) {
1399566063dSJacob Faibussowitsch     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
1409566063dSJacob Faibussowitsch     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
141c4762a1bSJed Brown   }
142c4762a1bSJed Brown   PetscFunctionReturn(0);
143c4762a1bSJed Brown }
144c4762a1bSJed Brown 
145c4762a1bSJed Brown /*
146c4762a1bSJed Brown      Defines the exact (analytic) solution to the ODE
147c4762a1bSJed Brown */
148*9371c9d4SSatish Balay static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx) {
149c4762a1bSJed Brown   const PetscScalar *uinit;
150c4762a1bSJed Brown   PetscScalar       *u, d0, q;
151c4762a1bSJed Brown 
152c4762a1bSJed Brown   PetscFunctionBegin;
1539566063dSJacob Faibussowitsch   PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit));
1549566063dSJacob Faibussowitsch   PetscCall(VecGetArrayWrite(U, &u));
155c4762a1bSJed Brown   d0 = uinit[0] - uinit[1];
156c4762a1bSJed Brown   if (d0 == 0.0) q = ctx->k * t;
157c4762a1bSJed Brown   else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0;
158c4762a1bSJed Brown   u[0] = uinit[0] / (1.0 + uinit[1] * q);
159c4762a1bSJed Brown   u[1] = u[0] - d0;
160c4762a1bSJed Brown   u[2] = uinit[1] + uinit[2] - u[1];
1619566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayWrite(U, &u));
1629566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit));
163c4762a1bSJed Brown   PetscFunctionReturn(0);
164c4762a1bSJed Brown }
165c4762a1bSJed Brown 
166*9371c9d4SSatish Balay int main(int argc, char **argv) {
167c4762a1bSJed Brown   TS                ts; /* ODE integrator */
168c4762a1bSJed Brown   Vec               U;  /* solution will be stored here */
169c4762a1bSJed Brown   Mat               A;  /* Jacobian matrix */
170c4762a1bSJed Brown   PetscMPIInt       size;
171c4762a1bSJed Brown   PetscInt          n = 3;
172c4762a1bSJed Brown   AppCtx            ctx;
173c4762a1bSJed Brown   PetscScalar      *u;
174c4762a1bSJed Brown   const char *const names[] = {"U1", "U2", "U3", NULL};
175c4762a1bSJed Brown 
176c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
177c4762a1bSJed Brown      Initialize program
178c4762a1bSJed Brown      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
179327415f7SBarry Smith   PetscFunctionBeginUser;
1809566063dSJacob Faibussowitsch   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
1819566063dSJacob Faibussowitsch   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
1823c633725SBarry Smith   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
183c4762a1bSJed Brown 
184c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185c4762a1bSJed Brown     Create necessary matrix and vectors
186c4762a1bSJed Brown     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1879566063dSJacob Faibussowitsch   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
1889566063dSJacob Faibussowitsch   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
1899566063dSJacob Faibussowitsch   PetscCall(MatSetFromOptions(A));
1909566063dSJacob Faibussowitsch   PetscCall(MatSetUp(A));
191c4762a1bSJed Brown 
1929566063dSJacob Faibussowitsch   PetscCall(MatCreateVecs(A, &U, NULL));
193c4762a1bSJed Brown 
194c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
195c4762a1bSJed Brown     Set runtime options
196c4762a1bSJed Brown     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
197c4762a1bSJed Brown   ctx.k = .9;
1989566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL));
1999566063dSJacob Faibussowitsch   PetscCall(VecDuplicate(U, &ctx.initialsolution));
2009566063dSJacob Faibussowitsch   PetscCall(VecGetArrayWrite(ctx.initialsolution, &u));
201c4762a1bSJed Brown   u[0] = 1;
202c4762a1bSJed Brown   u[1] = .7;
203c4762a1bSJed Brown   u[2] = 0;
2049566063dSJacob Faibussowitsch   PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u));
2059566063dSJacob Faibussowitsch   PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL));
206c4762a1bSJed Brown 
207c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
208c4762a1bSJed Brown      Create timestepping solver context
209c4762a1bSJed Brown      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2109566063dSJacob Faibussowitsch   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
2119566063dSJacob Faibussowitsch   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
2129566063dSJacob Faibussowitsch   PetscCall(TSSetType(ts, TSROSW));
2139566063dSJacob Faibussowitsch   PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &ctx));
2149566063dSJacob Faibussowitsch   PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobian)IJacobian, &ctx));
2159566063dSJacob Faibussowitsch   PetscCall(TSSetSolutionFunction(ts, (TSSolutionFunction)Solution, &ctx));
216c4762a1bSJed Brown 
217c4762a1bSJed Brown   {
218c4762a1bSJed Brown     DM    dm;
219c4762a1bSJed Brown     void *ptr;
2209566063dSJacob Faibussowitsch     PetscCall(TSGetDM(ts, &dm));
2219566063dSJacob Faibussowitsch     PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr));
2229566063dSJacob Faibussowitsch     PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr));
2239566063dSJacob Faibussowitsch     PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad));
2249566063dSJacob Faibussowitsch     PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad));
225c4762a1bSJed Brown   }
226c4762a1bSJed Brown 
227c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
228c4762a1bSJed Brown      Set initial conditions
229c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2309566063dSJacob Faibussowitsch   PetscCall(Solution(ts, 0, U, &ctx));
2319566063dSJacob Faibussowitsch   PetscCall(TSSetSolution(ts, U));
232c4762a1bSJed Brown 
233c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
234c4762a1bSJed Brown      Set solver options
235c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2369566063dSJacob Faibussowitsch   PetscCall(TSSetTimeStep(ts, .001));
2379566063dSJacob Faibussowitsch   PetscCall(TSSetMaxSteps(ts, 1000));
2389566063dSJacob Faibussowitsch   PetscCall(TSSetMaxTime(ts, 20.0));
2399566063dSJacob Faibussowitsch   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
2409566063dSJacob Faibussowitsch   PetscCall(TSSetFromOptions(ts));
2419566063dSJacob Faibussowitsch   PetscCall(TSMonitorLGSetVariableNames(ts, names));
242c4762a1bSJed Brown 
243c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244c4762a1bSJed Brown      Solve nonlinear system
245c4762a1bSJed Brown      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2469566063dSJacob Faibussowitsch   PetscCall(TSSolve(ts, U));
247c4762a1bSJed Brown 
2489566063dSJacob Faibussowitsch   PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD));
249c4762a1bSJed Brown 
250c4762a1bSJed Brown   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251c4762a1bSJed Brown      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
252c4762a1bSJed Brown    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
2539566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&ctx.initialsolution));
2549566063dSJacob Faibussowitsch   PetscCall(MatDestroy(&A));
2559566063dSJacob Faibussowitsch   PetscCall(VecDestroy(&U));
2569566063dSJacob Faibussowitsch   PetscCall(TSDestroy(&ts));
257c4762a1bSJed Brown 
2589566063dSJacob Faibussowitsch   PetscCall(PetscFinalize());
259b122ec5aSJacob Faibussowitsch   return 0;
260c4762a1bSJed Brown }
261c4762a1bSJed Brown 
262c4762a1bSJed Brown /*TEST
263c4762a1bSJed Brown 
264c4762a1bSJed Brown    test:
265c4762a1bSJed Brown      args: -ts_view
266dfd57a17SPierre Jolivet      requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
267c4762a1bSJed Brown 
268c4762a1bSJed Brown    test:
269c4762a1bSJed Brown      suffix: 2
270c4762a1bSJed Brown      args: -ts_monitor_lg_error -ts_monitor_lg_solution  -ts_view
271dfd57a17SPierre Jolivet      requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES)
272c4762a1bSJed Brown      output_file: output/ex1_1.out
273c4762a1bSJed Brown 
274c4762a1bSJed Brown TEST*/
275