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 PetscFunctionBegin; 80 PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR)); 81 PetscFunctionReturn(0); 82 } 83 84 PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v) { 85 PetscFunctionBegin; 86 PetscCall(PetscNew(ctx)); 87 PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR)); 88 PetscFunctionReturn(0); 89 } 90 91 /* 92 Defines the ODE passed to the ODE solver 93 */ 94 PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx) { 95 PetscScalar *f; 96 const PetscScalar *u, *udot; 97 98 PetscFunctionBegin; 99 /* The next three lines allow us to access the entries of the vectors directly */ 100 PetscCall(VecGetArrayRead(U, &u)); 101 PetscCall(VecGetArrayRead(Udot, &udot)); 102 PetscCall(VecGetArrayWrite(F, &f)); 103 f[0] = udot[0] + ctx->k * u[0] * u[1]; 104 f[1] = udot[1] + ctx->k * u[0] * u[1]; 105 f[2] = udot[2] - ctx->k * u[0] * u[1]; 106 PetscCall(VecRestoreArrayRead(U, &u)); 107 PetscCall(VecRestoreArrayRead(Udot, &udot)); 108 PetscCall(VecRestoreArrayWrite(F, &f)); 109 PetscFunctionReturn(0); 110 } 111 112 /* 113 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 114 */ 115 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx) { 116 PetscInt rowcol[] = {0, 1, 2}; 117 PetscScalar J[3][3]; 118 const PetscScalar *u, *udot; 119 120 PetscFunctionBegin; 121 PetscCall(VecGetArrayRead(U, &u)); 122 PetscCall(VecGetArrayRead(Udot, &udot)); 123 J[0][0] = a + ctx->k * u[1]; 124 J[0][1] = ctx->k * u[0]; 125 J[0][2] = 0.0; 126 J[1][0] = ctx->k * u[1]; 127 J[1][1] = a + ctx->k * u[0]; 128 J[1][2] = 0.0; 129 J[2][0] = -ctx->k * u[1]; 130 J[2][1] = -ctx->k * u[0]; 131 J[2][2] = a; 132 PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES)); 133 PetscCall(VecRestoreArrayRead(U, &u)); 134 PetscCall(VecRestoreArrayRead(Udot, &udot)); 135 136 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 137 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 138 if (A != B) { 139 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 140 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 141 } 142 PetscFunctionReturn(0); 143 } 144 145 /* 146 Defines the exact (analytic) solution to the ODE 147 */ 148 static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx) { 149 const PetscScalar *uinit; 150 PetscScalar *u, d0, q; 151 152 PetscFunctionBegin; 153 PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit)); 154 PetscCall(VecGetArrayWrite(U, &u)); 155 d0 = uinit[0] - uinit[1]; 156 if (d0 == 0.0) q = ctx->k * t; 157 else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0; 158 u[0] = uinit[0] / (1.0 + uinit[1] * q); 159 u[1] = u[0] - d0; 160 u[2] = uinit[1] + uinit[2] - u[1]; 161 PetscCall(VecRestoreArrayWrite(U, &u)); 162 PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit)); 163 PetscFunctionReturn(0); 164 } 165 166 int main(int argc, char **argv) { 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