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 PetscErrorCode ierr; 81 82 PetscFunctionBegin; 83 ierr = PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR);CHKERRQ(ierr); 84 PetscFunctionReturn(0); 85 } 86 87 PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v) 88 { 89 PetscErrorCode ierr; 90 91 PetscFunctionBegin; 92 ierr = PetscNew(ctx);CHKERRQ(ierr); 93 ierr = PetscViewerBinaryRead(v,&(*ctx)->k,1,NULL,PETSC_SCALAR);CHKERRQ(ierr); 94 PetscFunctionReturn(0); 95 } 96 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; 105 106 PetscFunctionBegin; 107 /* The next three lines allow us to access the entries of the vectors directly */ 108 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 109 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 110 ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 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 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 115 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 116 ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 117 PetscFunctionReturn(0); 118 } 119 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; 129 130 PetscFunctionBegin; 131 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 132 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 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 ierr = MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 137 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 138 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 139 140 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 141 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 142 if (A != B) { 143 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 144 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 145 } 146 PetscFunctionReturn(0); 147 } 148 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; 157 158 PetscFunctionBegin; 159 ierr = VecGetArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 160 ierr = VecGetArrayWrite(U,&u);CHKERRQ(ierr); 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 ierr = VecRestoreArrayWrite(U,&u);CHKERRQ(ierr); 168 ierr = VecRestoreArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 169 PetscFunctionReturn(0); 170 } 171 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 */ 177 PetscErrorCode ierr; 178 PetscMPIInt size; 179 PetscInt n = 3; 180 AppCtx ctx; 181 PetscScalar *u; 182 const char * const names[] = {"U1","U2","U3",NULL}; 183 184 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185 Initialize program 186 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 188 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 189 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 190 191 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 192 Create necessary matrix and vectors 193 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 194 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 195 ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 196 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 197 ierr = MatSetUp(A);CHKERRQ(ierr); 198 199 ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 200 201 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202 Set runtime options 203 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 204 ctx.k = .9; 205 ierr = PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL);CHKERRQ(ierr); 206 ierr = VecDuplicate(U,&ctx.initialsolution);CHKERRQ(ierr); 207 ierr = VecGetArrayWrite(ctx.initialsolution,&u);CHKERRQ(ierr); 208 u[0] = 1; 209 u[1] = .7; 210 u[2] = 0; 211 ierr = VecRestoreArrayWrite(ctx.initialsolution,&u);CHKERRQ(ierr); 212 ierr = PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL);CHKERRQ(ierr); 213 214 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215 Create timestepping solver context 216 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 218 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 219 ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); 220 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx);CHKERRQ(ierr); 221 ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);CHKERRQ(ierr); 222 ierr = TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx);CHKERRQ(ierr); 223 224 { 225 DM dm; 226 void *ptr; 227 ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 228 ierr = PetscDLSym(NULL,"IFunctionView",&ptr);CHKERRQ(ierr); 229 ierr = PetscDLSym(NULL,"IFunctionLoad",&ptr);CHKERRQ(ierr); 230 ierr = DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 231 ierr = DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 232 } 233 234 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 235 Set initial conditions 236 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 237 ierr = Solution(ts,0,U,&ctx);CHKERRQ(ierr); 238 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 239 240 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 241 Set solver options 242 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 243 ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr); 244 ierr = TSSetMaxSteps(ts,1000);CHKERRQ(ierr); 245 ierr = TSSetMaxTime(ts,20.0);CHKERRQ(ierr); 246 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 247 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 248 ierr = TSMonitorLGSetVariableNames(ts,names);CHKERRQ(ierr); 249 250 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 251 Solve nonlinear system 252 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 253 ierr = TSSolve(ts,U);CHKERRQ(ierr); 254 255 ierr = TSView(ts,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 256 257 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 258 Free work space. All PETSc objects should be destroyed when they are no longer needed. 259 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 260 ierr = VecDestroy(&ctx.initialsolution);CHKERRQ(ierr); 261 ierr = MatDestroy(&A);CHKERRQ(ierr); 262 ierr = VecDestroy(&U);CHKERRQ(ierr); 263 ierr = TSDestroy(&ts);CHKERRQ(ierr); 264 265 ierr = PetscFinalize(); 266 return ierr; 267 } 268 269 270 /*TEST 271 272 test: 273 args: -ts_view 274 requires: dlsym define(PETSC_HAVE_DYNAMIC_LIBRARIES) 275 276 test: 277 suffix: 2 278 args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 279 requires: x dlsym define(PETSC_HAVE_DYNAMIC_LIBRARIES) 280 output_file: output/ex1_1.out 281 282 TEST*/ 283