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 PetscErrorCode ierr; 171 PetscMPIInt size; 172 PetscInt n = 3; 173 AppCtx ctx; 174 PetscScalar *u; 175 const char * const names[] = {"U1","U2","U3",NULL}; 176 177 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178 Initialize program 179 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 180 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 181 CHKERRMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 182 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 183 184 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185 Create necessary matrix and vectors 186 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187 CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 188 CHKERRQ(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 189 CHKERRQ(MatSetFromOptions(A)); 190 CHKERRQ(MatSetUp(A)); 191 192 CHKERRQ(MatCreateVecs(A,&U,NULL)); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Set runtime options 196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197 ctx.k = .9; 198 CHKERRQ(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 199 CHKERRQ(VecDuplicate(U,&ctx.initialsolution)); 200 CHKERRQ(VecGetArrayWrite(ctx.initialsolution,&u)); 201 u[0] = 1; 202 u[1] = .7; 203 u[2] = 0; 204 CHKERRQ(VecRestoreArrayWrite(ctx.initialsolution,&u)); 205 CHKERRQ(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 206 207 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 208 Create timestepping solver context 209 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 210 CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 211 CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 212 CHKERRQ(TSSetType(ts,TSROSW)); 213 CHKERRQ(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 214 CHKERRQ(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 215 CHKERRQ(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 216 217 { 218 DM dm; 219 void *ptr; 220 CHKERRQ(TSGetDM(ts,&dm)); 221 CHKERRQ(PetscDLSym(NULL,"IFunctionView",&ptr)); 222 CHKERRQ(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 223 CHKERRQ(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 224 CHKERRQ(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 225 } 226 227 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 228 Set initial conditions 229 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 230 CHKERRQ(Solution(ts,0,U,&ctx)); 231 CHKERRQ(TSSetSolution(ts,U)); 232 233 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 234 Set solver options 235 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 236 CHKERRQ(TSSetTimeStep(ts,.001)); 237 CHKERRQ(TSSetMaxSteps(ts,1000)); 238 CHKERRQ(TSSetMaxTime(ts,20.0)); 239 CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 240 CHKERRQ(TSSetFromOptions(ts)); 241 CHKERRQ(TSMonitorLGSetVariableNames(ts,names)); 242 243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244 Solve nonlinear system 245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 246 CHKERRQ(TSSolve(ts,U)); 247 248 CHKERRQ(TSView(ts,PETSC_VIEWER_BINARY_WORLD)); 249 250 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 251 Free work space. All PETSc objects should be destroyed when they are no longer needed. 252 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 253 CHKERRQ(VecDestroy(&ctx.initialsolution)); 254 CHKERRQ(MatDestroy(&A)); 255 CHKERRQ(VecDestroy(&U)); 256 CHKERRQ(TSDestroy(&ts)); 257 258 ierr = PetscFinalize(); 259 return ierr; 260 } 261 262 /*TEST 263 264 test: 265 args: -ts_view 266 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 267 268 test: 269 suffix: 2 270 args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 271 requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 272 output_file: output/ex1_1.out 273 274 TEST*/ 275