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 progress 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 formatted 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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(PETSC_SUCCESS); 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]; 128 J[0][1] = ctx->k * u[0]; 129 J[0][2] = 0.0; 130 J[1][0] = ctx->k * u[1]; 131 J[1][1] = a + ctx->k * u[0]; 132 J[1][2] = 0.0; 133 J[2][0] = -ctx->k * u[1]; 134 J[2][1] = -ctx->k * u[0]; 135 J[2][2] = a; 136 PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES)); 137 PetscCall(VecRestoreArrayRead(U, &u)); 138 PetscCall(VecRestoreArrayRead(Udot, &udot)); 139 140 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 141 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 142 if (A != B) { 143 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 144 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 145 } 146 PetscFunctionReturn(PETSC_SUCCESS); 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 const PetscScalar *uinit; 155 PetscScalar *u, d0, q; 156 157 PetscFunctionBegin; 158 PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit)); 159 PetscCall(VecGetArrayWrite(U, &u)); 160 d0 = uinit[0] - uinit[1]; 161 if (d0 == 0.0) q = ctx->k * t; 162 else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0; 163 u[0] = uinit[0] / (1.0 + uinit[1] * q); 164 u[1] = u[0] - d0; 165 u[2] = uinit[1] + uinit[2] - u[1]; 166 PetscCall(VecRestoreArrayWrite(U, &u)); 167 PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit)); 168 PetscFunctionReturn(PETSC_SUCCESS); 169 } 170 171 int main(int argc, char **argv) 172 { 173 TS ts; /* ODE integrator */ 174 Vec U; /* solution will be stored here */ 175 Mat A; /* Jacobian matrix */ 176 PetscMPIInt size; 177 PetscInt n = 3; 178 AppCtx ctx; 179 PetscScalar *u; 180 const char *const names[] = {"U1", "U2", "U3", NULL}; 181 182 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 183 Initialize program 184 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 185 PetscFunctionBeginUser; 186 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 187 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 188 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 189 190 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 191 Create necessary matrix and vectors 192 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 193 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 194 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE)); 195 PetscCall(MatSetFromOptions(A)); 196 PetscCall(MatSetUp(A)); 197 198 PetscCall(MatCreateVecs(A, &U, NULL)); 199 200 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 201 Set runtime options 202 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 203 ctx.k = .9; 204 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL)); 205 PetscCall(VecDuplicate(U, &ctx.initialsolution)); 206 PetscCall(VecGetArrayWrite(ctx.initialsolution, &u)); 207 u[0] = 1; 208 u[1] = .7; 209 u[2] = 0; 210 PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u)); 211 PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL)); 212 213 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 214 Create timestepping solver context 215 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 216 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 217 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 218 PetscCall(TSSetType(ts, TSROSW)); 219 PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &ctx)); 220 PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobian)IJacobian, &ctx)); 221 PetscCall(TSSetSolutionFunction(ts, (TSSolutionFunction)Solution, &ctx)); 222 223 { 224 DM dm; 225 void *ptr; 226 PetscCall(TSGetDM(ts, &dm)); 227 PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr)); 228 PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr)); 229 PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 230 PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 231 } 232 233 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 234 Set initial conditions 235 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 236 PetscCall(Solution(ts, 0, U, &ctx)); 237 PetscCall(TSSetSolution(ts, U)); 238 239 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 240 Set solver options 241 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 242 PetscCall(TSSetTimeStep(ts, .001)); 243 PetscCall(TSSetMaxSteps(ts, 1000)); 244 PetscCall(TSSetMaxTime(ts, 20.0)); 245 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 246 PetscCall(TSSetFromOptions(ts)); 247 PetscCall(TSMonitorLGSetVariableNames(ts, names)); 248 249 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 250 Solve nonlinear system 251 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 252 PetscCall(TSSolve(ts, U)); 253 254 PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD)); 255 256 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 257 Free work space. All PETSc objects should be destroyed when they are no longer needed. 258 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 259 PetscCall(VecDestroy(&ctx.initialsolution)); 260 PetscCall(MatDestroy(&A)); 261 PetscCall(VecDestroy(&U)); 262 PetscCall(TSDestroy(&ts)); 263 264 PetscCall(PetscFinalize()); 265 return 0; 266 } 267 268 /*TEST 269 270 test: 271 args: -ts_view 272 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 273 274 test: 275 suffix: 2 276 args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 277 requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 278 output_file: output/ex1_1.out 279 280 TEST*/ 281