1 2 static char help[] = "Nonlinear Reaction Problem from Chemistry.\n"; 3 4 /*F 5 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 9 10 Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry 11 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} 17 18 Helpful runtime monitoring options: 19 -ts_view - prints information about the solver being used 20 -ts_monitor - prints the progess 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 26 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 32 33 F*/ 34 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> 41 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 44 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> 58 59 */ 60 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 */ 70 71 #include <petscts.h> 72 73 typedef struct { 74 PetscScalar k; 75 Vec initialsolution; 76 } AppCtx; 77 78 PetscErrorCode IFunctionView(AppCtx *ctx,PetscViewer v) 79 { 80 PetscFunctionBegin; 81 CHKERRQ(PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR)); 82 PetscFunctionReturn(0); 83 } 84 85 PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v) 86 { 87 PetscFunctionBegin; 88 CHKERRQ(PetscNew(ctx)); 89 CHKERRQ(PetscViewerBinaryRead(v,&(*ctx)->k,1,NULL,PETSC_SCALAR)); 90 PetscFunctionReturn(0); 91 } 92 93 /* 94 Defines the ODE passed to the ODE solver 95 */ 96 PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 97 { 98 PetscScalar *f; 99 const PetscScalar *u,*udot; 100 101 PetscFunctionBegin; 102 /* The next three lines allow us to access the entries of the vectors directly */ 103 CHKERRQ(VecGetArrayRead(U,&u)); 104 CHKERRQ(VecGetArrayRead(Udot,&udot)); 105 CHKERRQ(VecGetArrayWrite(F,&f)); 106 f[0] = udot[0] + ctx->k*u[0]*u[1]; 107 f[1] = udot[1] + ctx->k*u[0]*u[1]; 108 f[2] = udot[2] - ctx->k*u[0]*u[1]; 109 CHKERRQ(VecRestoreArrayRead(U,&u)); 110 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 111 CHKERRQ(VecRestoreArrayWrite(F,&f)); 112 PetscFunctionReturn(0); 113 } 114 115 /* 116 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 117 */ 118 PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx) 119 { 120 PetscInt rowcol[] = {0,1,2}; 121 PetscScalar J[3][3]; 122 const PetscScalar *u,*udot; 123 124 PetscFunctionBegin; 125 CHKERRQ(VecGetArrayRead(U,&u)); 126 CHKERRQ(VecGetArrayRead(Udot,&udot)); 127 J[0][0] = a + ctx->k*u[1]; J[0][1] = ctx->k*u[0]; J[0][2] = 0.0; 128 J[1][0] = ctx->k*u[1]; J[1][1] = a + ctx->k*u[0]; J[1][2] = 0.0; 129 J[2][0] = -ctx->k*u[1]; J[2][1] = -ctx->k*u[0]; J[2][2] = a; 130 CHKERRQ(MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES)); 131 CHKERRQ(VecRestoreArrayRead(U,&u)); 132 CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 133 134 CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 135 CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 136 if (A != B) { 137 CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 138 CHKERRQ(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 139 } 140 PetscFunctionReturn(0); 141 } 142 143 /* 144 Defines the exact (analytic) solution to the ODE 145 */ 146 static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx) 147 { 148 const PetscScalar *uinit; 149 PetscScalar *u,d0,q; 150 151 PetscFunctionBegin; 152 CHKERRQ(VecGetArrayRead(ctx->initialsolution,&uinit)); 153 CHKERRQ(VecGetArrayWrite(U,&u)); 154 d0 = uinit[0] - uinit[1]; 155 if (d0 == 0.0) q = ctx->k*t; 156 else q = (1.0 - PetscExpScalar(-ctx->k*t*d0))/d0; 157 u[0] = uinit[0]/(1.0 + uinit[1]*q); 158 u[1] = u[0] - d0; 159 u[2] = uinit[1] + uinit[2] - u[1]; 160 CHKERRQ(VecRestoreArrayWrite(U,&u)); 161 CHKERRQ(VecRestoreArrayRead(ctx->initialsolution,&uinit)); 162 PetscFunctionReturn(0); 163 } 164 165 int main(int argc,char **argv) 166 { 167 TS ts; /* ODE integrator */ 168 Vec U; /* solution will be stored here */ 169 Mat A; /* Jacobian matrix */ 170 PetscMPIInt size; 171 PetscInt n = 3; 172 AppCtx ctx; 173 PetscScalar *u; 174 const char * const names[] = {"U1","U2","U3",NULL}; 175 176 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 177 Initialize program 178 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 179 CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 180 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 181 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 182 183 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184 Create necessary matrix and vectors 185 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 186 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 187 CHKERRQ(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 188 CHKERRQ(MatSetFromOptions(A)); 189 CHKERRQ(MatSetUp(A)); 190 191 CHKERRQ(MatCreateVecs(A,&U,NULL)); 192 193 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194 Set runtime options 195 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 196 ctx.k = .9; 197 CHKERRQ(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 198 CHKERRQ(VecDuplicate(U,&ctx.initialsolution)); 199 CHKERRQ(VecGetArrayWrite(ctx.initialsolution,&u)); 200 u[0] = 1; 201 u[1] = .7; 202 u[2] = 0; 203 CHKERRQ(VecRestoreArrayWrite(ctx.initialsolution,&u)); 204 CHKERRQ(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 205 206 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207 Create timestepping solver context 208 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 209 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 210 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 211 CHKERRQ(TSSetType(ts,TSROSW)); 212 CHKERRQ(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 213 CHKERRQ(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 214 CHKERRQ(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 215 216 { 217 DM dm; 218 void *ptr; 219 CHKERRQ(TSGetDM(ts,&dm)); 220 CHKERRQ(PetscDLSym(NULL,"IFunctionView",&ptr)); 221 CHKERRQ(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 222 CHKERRQ(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 223 CHKERRQ(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 224 } 225 226 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 227 Set initial conditions 228 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 229 CHKERRQ(Solution(ts,0,U,&ctx)); 230 CHKERRQ(TSSetSolution(ts,U)); 231 232 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233 Set solver options 234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 235 CHKERRQ(TSSetTimeStep(ts,.001)); 236 CHKERRQ(TSSetMaxSteps(ts,1000)); 237 CHKERRQ(TSSetMaxTime(ts,20.0)); 238 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 239 CHKERRQ(TSSetFromOptions(ts)); 240 CHKERRQ(TSMonitorLGSetVariableNames(ts,names)); 241 242 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243 Solve nonlinear system 244 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 245 CHKERRQ(TSSolve(ts,U)); 246 247 CHKERRQ(TSView(ts,PETSC_VIEWER_BINARY_WORLD)); 248 249 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 250 Free work space. All PETSc objects should be destroyed when they are no longer needed. 251 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 252 CHKERRQ(VecDestroy(&ctx.initialsolution)); 253 CHKERRQ(MatDestroy(&A)); 254 CHKERRQ(VecDestroy(&U)); 255 CHKERRQ(TSDestroy(&ts)); 256 257 CHKERRQ(PetscFinalize()); 258 return 0; 259 } 260 261 /*TEST 262 263 test: 264 args: -ts_view 265 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 266 267 test: 268 suffix: 2 269 args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 270 requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 271 output_file: output/ex1_1.out 272 273 TEST*/ 274