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*/