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_$_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 PetscFunctionBeginUser; 180 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 181 PetscCallMPI(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 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 188 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 189 PetscCall(MatSetFromOptions(A)); 190 PetscCall(MatSetUp(A)); 191 192 PetscCall(MatCreateVecs(A,&U,NULL)); 193 194 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195 Set runtime options 196 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197 ctx.k = .9; 198 PetscCall(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 199 PetscCall(VecDuplicate(U,&ctx.initialsolution)); 200 PetscCall(VecGetArrayWrite(ctx.initialsolution,&u)); 201 u[0] = 1; 202 u[1] = .7; 203 u[2] = 0; 204 PetscCall(VecRestoreArrayWrite(ctx.initialsolution,&u)); 205 PetscCall(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 206 207 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 208 Create timestepping solver context 209 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 210 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 211 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 212 PetscCall(TSSetType(ts,TSROSW)); 213 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 214 PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 215 PetscCall(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 216 217 { 218 DM dm; 219 void *ptr; 220 PetscCall(TSGetDM(ts,&dm)); 221 PetscCall(PetscDLSym(NULL,"IFunctionView",&ptr)); 222 PetscCall(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 223 PetscCall(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 224 PetscCall(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 225 } 226 227 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 228 Set initial conditions 229 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 230 PetscCall(Solution(ts,0,U,&ctx)); 231 PetscCall(TSSetSolution(ts,U)); 232 233 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 234 Set solver options 235 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 236 PetscCall(TSSetTimeStep(ts,.001)); 237 PetscCall(TSSetMaxSteps(ts,1000)); 238 PetscCall(TSSetMaxTime(ts,20.0)); 239 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 240 PetscCall(TSSetFromOptions(ts)); 241 PetscCall(TSMonitorLGSetVariableNames(ts,names)); 242 243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244 Solve nonlinear system 245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 246 PetscCall(TSSolve(ts,U)); 247 248 PetscCall(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 PetscCall(VecDestroy(&ctx.initialsolution)); 254 PetscCall(MatDestroy(&A)); 255 PetscCall(VecDestroy(&U)); 256 PetscCall(TSDestroy(&ts)); 257 258 PetscCall(PetscFinalize()); 259 return 0; 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