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 PetscCall(PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR)); 82 PetscFunctionReturn(0); 83 } 84 85 PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v) 86 { 87 PetscFunctionBegin; 88 PetscCall(PetscNew(ctx)); 89 PetscCall(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 PetscCall(VecGetArrayRead(U,&u)); 104 PetscCall(VecGetArrayRead(Udot,&udot)); 105 PetscCall(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 PetscCall(VecRestoreArrayRead(U,&u)); 110 PetscCall(VecRestoreArrayRead(Udot,&udot)); 111 PetscCall(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 PetscCall(VecGetArrayRead(U,&u)); 126 PetscCall(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 PetscCall(MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES)); 131 PetscCall(VecRestoreArrayRead(U,&u)); 132 PetscCall(VecRestoreArrayRead(Udot,&udot)); 133 134 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 135 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 136 if (A != B) { 137 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 138 PetscCall(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 PetscCall(VecGetArrayRead(ctx->initialsolution,&uinit)); 153 PetscCall(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 PetscCall(VecRestoreArrayWrite(U,&u)); 161 PetscCall(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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 180 PetscCallMPI(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 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 187 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 188 PetscCall(MatSetFromOptions(A)); 189 PetscCall(MatSetUp(A)); 190 191 PetscCall(MatCreateVecs(A,&U,NULL)); 192 193 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194 Set runtime options 195 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 196 ctx.k = .9; 197 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 198 PetscCall(VecDuplicate(U,&ctx.initialsolution)); 199 PetscCall(VecGetArrayWrite(ctx.initialsolution,&u)); 200 u[0] = 1; 201 u[1] = .7; 202 u[2] = 0; 203 PetscCall(VecRestoreArrayWrite(ctx.initialsolution,&u)); 204 PetscCall(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 205 206 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207 Create timestepping solver context 208 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 209 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 210 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 211 PetscCall(TSSetType(ts,TSROSW)); 212 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 213 PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 214 PetscCall(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 215 216 { 217 DM dm; 218 void *ptr; 219 PetscCall(TSGetDM(ts,&dm)); 220 PetscCall(PetscDLSym(NULL,"IFunctionView",&ptr)); 221 PetscCall(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 222 PetscCall(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 223 PetscCall(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 224 } 225 226 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 227 Set initial conditions 228 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 229 PetscCall(Solution(ts,0,U,&ctx)); 230 PetscCall(TSSetSolution(ts,U)); 231 232 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233 Set solver options 234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 235 PetscCall(TSSetTimeStep(ts,.001)); 236 PetscCall(TSSetMaxSteps(ts,1000)); 237 PetscCall(TSSetMaxTime(ts,20.0)); 238 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 239 PetscCall(TSSetFromOptions(ts)); 240 PetscCall(TSMonitorLGSetVariableNames(ts,names)); 241 242 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243 Solve nonlinear system 244 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 245 PetscCall(TSSolve(ts,U)); 246 247 PetscCall(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 PetscCall(VecDestroy(&ctx.initialsolution)); 253 PetscCall(MatDestroy(&A)); 254 PetscCall(VecDestroy(&U)); 255 PetscCall(TSDestroy(&ts)); 256 257 PetscCall(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