Actual source code: ex1.c
petsc-3.11.1 2019-04-12
2: static char help[] = " Nonlinear Reaction Problem from Chemistry.\n" ;
This directory contains examples based on the PDES/ODES given in the book
Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by
W. Hundsdorf and J.G. Verwer
Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry
\begin{eqnarray}
{U_1}_t - k U_1 U_2 & = & 0 \\
{U_2}_t - k U_1 U_2 & = & 0 \\
{U_3}_t - k U_1 U_2 & = & 0
\end{eqnarray}
Helpful runtime monitoring options:
-ts_view - prints information about the solver being used
-ts_monitor - prints the progess of the solver
-ts_adapt_monitor - prints the progress of the time-step adaptor
-ts_monitor_lg_timestep - plots the size of each timestep ( at each time-step)
-ts_monitor_lg_solution - plots each component of the solution as a function of time ( at each timestep)
-ts_monitor_lg_error - plots each component of the error in the solution as a function of time ( at each timestep)
-draw_pause -2 - hold the plots a the end of the solution process, enter a mouse press in each window to end the process
-ts_monitor_lg_timestep -1 - plots the size of each timestep ( at the end of the solution process)
-ts_monitor_lg_solution -1 - plots each component of the solution as a function of time ( at the end of the solution process)
-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)
-lg_use_markers false - do NOT show the data points on the plots
-draw_save - save the timestep and solution plot as a .Gif image file
35: /*
36: Project: Generate a nicely formated 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_ZZZ_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: {
83: PetscViewerBinaryWrite ( v,&ctx->k,1,PETSC_SCALAR,PETSC_FALSE );
84: return ( 0);
85: }
87: PetscErrorCode IFunctionLoad( AppCtx **ctx,PetscViewer v)
88: {
92: PetscNew ( ctx);
93: PetscViewerBinaryRead ( v,&( *ctx)->k,1,NULL,PETSC_SCALAR);
94: return ( 0);
95: }
97: /*
98: Defines the ODE passed to the ODE solver
99: */
100: PetscErrorCode IFunction( TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx)
101: {
102: PetscErrorCode ierr;
103: PetscScalar *f;
104: const PetscScalar *u,*udot;
107: /* The next three lines allow us to access the entries of the vectors directly */
108: VecGetArrayRead ( U,&u);
109: VecGetArrayRead ( Udot,&udot);
110: VecGetArray ( F,&f);
111: f[0] = udot[0] + ctx->k*u[0]*u[1];
112: f[1] = udot[1] + ctx->k*u[0]*u[1];
113: f[2] = udot[2] - ctx->k*u[0]*u[1];
114: VecRestoreArrayRead ( U,&u);
115: VecRestoreArrayRead ( Udot,&udot);
116: VecRestoreArray ( F,&f);
117: return ( 0);
118: }
120: /*
121: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian ( ) for the meaning of a and the Jacobian.
122: */
123: PetscErrorCode IJacobian( TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx)
124: {
125: PetscErrorCode ierr;
126: PetscInt rowcol[] = {0,1,2};
127: PetscScalar J[3][3];
128: const PetscScalar *u,*udot;
131: VecGetArrayRead ( U,&u);
132: VecGetArrayRead ( Udot,&udot);
133: J[0][0] = a + ctx->k*u[1]; J[0][1] = ctx->k*u[0]; J[0][2] = 0.0;
134: J[1][0] = ctx->k*u[1]; J[1][1] = a + ctx->k*u[0]; J[1][2] = 0.0;
135: J[2][0] = -ctx->k*u[1]; J[2][1] = -ctx->k*u[0]; J[2][2] = a;
136: MatSetValues ( B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES );
137: VecRestoreArrayRead ( U,&u);
138: VecRestoreArrayRead ( Udot,&udot);
140: MatAssemblyBegin ( A,MAT_FINAL_ASSEMBLY );
141: MatAssemblyEnd ( A,MAT_FINAL_ASSEMBLY );
142: if ( A != B) {
143: MatAssemblyBegin ( B,MAT_FINAL_ASSEMBLY );
144: MatAssemblyEnd ( B,MAT_FINAL_ASSEMBLY );
145: }
146: return ( 0);
147: }
149: /*
150: Defines the exact ( analytic) solution to the ODE
151: */
152: static PetscErrorCode Solution( TS ts,PetscReal t,Vec U,AppCtx *ctx)
153: {
154: PetscErrorCode ierr;
155: const PetscScalar *uinit;
156: PetscScalar *u,d0,q;
159: VecGetArrayRead ( ctx->initialsolution,&uinit);
160: VecGetArray ( U,&u);
161: d0 = uinit[0] - uinit[1];
162: if ( d0 = = 0.0) q = ctx->k*t;
163: else q = ( 1.0 - PetscExpScalar( -ctx->k*t*d0))/d0;
164: u[0] = uinit[0]/( 1.0 + uinit[1]*q);
165: u[1] = u[0] - d0;
166: u[2] = uinit[1] + uinit[2] - u[1];
167: VecRestoreArray ( U,&u);
168: VecRestoreArrayRead ( ctx->initialsolution,&uinit);
169: return ( 0);
170: }
172: int main( int argc,char **argv)
173: {
174: TS ts; /* ODE integrator */
175: Vec U; /* solution will be stored here */
176: Mat A; /* Jacobian matrix */
178: PetscMPIInt size;
179: PetscInt n = 3;
180: AppCtx ctx;
181: PetscScalar *u;
182: const char * const names[] = {" U1" ," U2" ," U3" ,NULL};
184: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
185: Initialize program
186: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
187: PetscInitialize ( &argc,&argv,( char*)0,help);if ( ierr) return ierr;
188: MPI_Comm_size ( PETSC_COMM_WORLD ,&size);
189: if ( size > 1) SETERRQ ( PETSC_COMM_WORLD ,PETSC_ERR_SUP," Only for sequential runs" );
191: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
192: Create necessary matrix and vectors
193: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
194: MatCreate ( PETSC_COMM_WORLD ,&A);
195: MatSetSizes ( A,n,n,PETSC_DETERMINE ,PETSC_DETERMINE );
196: MatSetFromOptions ( A);
197: MatSetUp ( A);
199: MatCreateVecs ( A,&U,NULL);
201: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
202: Set runtime options
203: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
204: ctx.k = .9;
205: PetscOptionsGetScalar ( NULL,NULL," -k" ,&ctx.k,NULL);
206: VecDuplicate ( U,&ctx.initialsolution);
207: VecGetArray ( ctx.initialsolution,&u);
208: u[0] = 1;
209: u[1] = .7;
210: u[2] = 0;
211: VecRestoreArray ( ctx.initialsolution,&u);
212: PetscOptionsGetVec( NULL,NULL," -initial" ,ctx.initialsolution,NULL);
214: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
215: Create timestepping solver context
216: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
217: TSCreate ( PETSC_COMM_WORLD ,&ts);
218: TSSetProblemType ( ts,TS_NONLINEAR );
219: TSSetType ( ts,TSROSW );
220: TSSetIFunction ( ts,NULL,( TSIFunction) IFunction,&ctx);
221: TSSetIJacobian ( ts,A,A,( TSIJacobian)IJacobian,&ctx);
222: TSSetSolutionFunction ( ts,( TSSolutionFunction)Solution,&ctx);
224: {
225: DM dm;
226: void *ptr;
227: TSGetDM ( ts,&dm);
228: PetscDLSym ( NULL," IFunctionView" ,&ptr);
229: PetscDLSym ( NULL," IFunctionLoad" ,&ptr);
230: DMTSSetIFunctionSerialize ( dm,( PetscErrorCode ( *)( void*,PetscViewer ))IFunctionView,( PetscErrorCode ( *)( void**,PetscViewer ))IFunctionLoad);
231: DMTSSetIJacobianSerialize ( dm,( PetscErrorCode ( *)( void*,PetscViewer ))IFunctionView,( PetscErrorCode ( *)( void**,PetscViewer ))IFunctionLoad);
232: }
234: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
235: Set initial conditions
236: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
237: Solution( ts,0,U,&ctx);
238: TSSetSolution ( ts,U);
240: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
241: Set solver options
242: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
243: TSSetTimeStep ( ts,.001);
244: TSSetMaxSteps ( ts,1000);
245: TSSetMaxTime ( ts,20.0);
246: TSSetExactFinalTime ( ts,TS_EXACTFINALTIME_STEPOVER );
247: TSSetFromOptions ( ts);
248: TSMonitorLGSetVariableNames ( ts,names);
250: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
251: Solve nonlinear system
252: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
253: TSSolve ( ts,U);
255: TSView ( ts,PETSC_VIEWER_BINARY_WORLD );
257: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258: Free work space. All PETSc objects should be destroyed when they are no longer needed.
259: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260: VecDestroy ( &ctx.initialsolution);
261: MatDestroy ( &A);
262: VecDestroy ( &U);
263: TSDestroy ( &ts);
265: PetscFinalize ( );
266: return ierr;
267: }
270: /*TEST
272: test:
273: args: -ts_view
274: requires: dlsym define( PETSC_HAVE_DYNAMIC_LIBRARIES)
276: test:
277: suffix: 2
278: args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view
279: requires: x
280: output_file: output/ex1_1.out
281: requires: dlsym define( PETSC_HAVE_DYNAMIC_LIBRARIES)
283: TEST*/