1 static char help[] = "Nonlinear Reaction Problem from Chemistry.\n"; 2 3 /*F 4 5 This directory contains examples based on the PDES/ODES given in the book 6 Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by 7 W. Hundsdorf and J.G. Verwer 8 9 Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry 10 11 \begin{eqnarray} 12 {U_1}_t - k U_1 U_2 & = & 0 \\ 13 {U_2}_t - k U_1 U_2 & = & 0 \\ 14 {U_3}_t - k U_1 U_2 & = & 0 15 \end{eqnarray} 16 17 Helpful runtime monitoring options: 18 -ts_view - prints information about the solver being used 19 -ts_monitor - prints the progress of the solver 20 -ts_adapt_monitor - prints the progress of the time-step adaptor 21 -ts_monitor_lg_timestep - plots the size of each timestep (at each time-step) 22 -ts_monitor_lg_solution - plots each component of the solution as a function of time (at each timestep) 23 -ts_monitor_lg_error - plots each component of the error in the solution as a function of time (at each timestep) 24 -draw_pause -2 - hold the plots a the end of the solution process, enter a mouse press in each window to end the process 25 26 -ts_monitor_lg_timestep -1 - plots the size of each timestep (at the end of the solution process) 27 -ts_monitor_lg_solution -1 - plots each component of the solution as a function of time (at the end of the solution process) 28 -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) 29 -lg_use_markers false - do NOT show the data points on the plots 30 -draw_save - save the timestep and solution plot as a .Gif image file 31 32 F*/ 33 34 /* 35 Project: Generate a nicely formatted HTML page using 36 1) the HTML version of the source code and text in this file, use make html to generate the file ex1.c.html 37 2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_$_1_0.Gif) and 38 3) the text output (output.txt) generated by running the following commands. 39 4) <iframe src="generated_topics.html" scrolling="no" frameborder="0" width=600 height=300></iframe> 40 41 rm -rf *.Gif 42 ./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 43 44 For example something like 45 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> 46 <html> 47 <head> 48 <meta http-equiv="content-type" content="text/html;charset=utf-8"> 49 <title>PETSc Example -- Nonlinear Reaction Problem from Chemistry</title> 50 </head> 51 <body> 52 <iframe src="ex1.c.html" scrolling="yes" frameborder="1" width=2000 height=400></iframe> 53 <img alt="" src="Draw_0x84000000_0_0.Gif"/><img alt="" src="Draw_0x84000001_0_0.Gif"/><img alt="" src="Draw_0x84000001_1_0.Gif"/> 54 <iframe src="output.txt" scrolling="yes" frameborder="1" width=2000 height=1000></iframe> 55 </body> 56 </html> 57 58 */ 59 60 /* 61 Include "petscts.h" so that we can use TS solvers. Note that this 62 file automatically includes: 63 petscsys.h - base PETSc routines petscvec.h - vectors 64 petscmat.h - matrices 65 petscis.h - index sets petscksp.h - Krylov subspace methods 66 petscviewer.h - viewers petscpc.h - preconditioners 67 petscksp.h - linear solvers 68 */ 69 70 #include <petscts.h> 71 72 typedef struct { 73 PetscScalar k; 74 Vec initialsolution; 75 } AppCtx; 76 77 PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v) 78 { 79 PetscFunctionBegin; 80 PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR)); 81 PetscFunctionReturn(PETSC_SUCCESS); 82 } 83 84 PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v) 85 { 86 PetscFunctionBegin; 87 PetscCall(PetscNew(ctx)); 88 PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR)); 89 PetscFunctionReturn(PETSC_SUCCESS); 90 } 91 92 /* 93 Defines the ODE passed to the ODE solver 94 */ 95 PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx) 96 { 97 PetscScalar *f; 98 const PetscScalar *u, *udot; 99 100 PetscFunctionBegin; 101 /* The next three lines allow us to access the entries of the vectors directly */ 102 PetscCall(VecGetArrayRead(U, &u)); 103 PetscCall(VecGetArrayRead(Udot, &udot)); 104 PetscCall(VecGetArrayWrite(F, &f)); 105 f[0] = udot[0] + ctx->k * u[0] * u[1]; 106 f[1] = udot[1] + ctx->k * u[0] * u[1]; 107 f[2] = udot[2] - ctx->k * u[0] * u[1]; 108 PetscCall(VecRestoreArrayRead(U, &u)); 109 PetscCall(VecRestoreArrayRead(Udot, &udot)); 110 PetscCall(VecRestoreArrayWrite(F, &f)); 111 PetscFunctionReturn(PETSC_SUCCESS); 112 } 113 114 /* 115 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 116 */ 117 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx) 118 { 119 PetscInt rowcol[] = {0, 1, 2}; 120 PetscScalar J[3][3]; 121 const PetscScalar *u, *udot; 122 123 PetscFunctionBegin; 124 PetscCall(VecGetArrayRead(U, &u)); 125 PetscCall(VecGetArrayRead(Udot, &udot)); 126 J[0][0] = a + ctx->k * u[1]; 127 J[0][1] = ctx->k * u[0]; 128 J[0][2] = 0.0; 129 J[1][0] = ctx->k * u[1]; 130 J[1][1] = a + ctx->k * u[0]; 131 J[1][2] = 0.0; 132 J[2][0] = -ctx->k * u[1]; 133 J[2][1] = -ctx->k * u[0]; 134 J[2][2] = a; 135 PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES)); 136 PetscCall(VecRestoreArrayRead(U, &u)); 137 PetscCall(VecRestoreArrayRead(Udot, &udot)); 138 139 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 140 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 141 if (A != B) { 142 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 143 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 144 } 145 PetscFunctionReturn(PETSC_SUCCESS); 146 } 147 148 /* 149 Defines the exact (analytic) solution to the ODE 150 */ 151 static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx) 152 { 153 const PetscScalar *uinit; 154 PetscScalar *u, d0, q; 155 156 PetscFunctionBegin; 157 PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit)); 158 PetscCall(VecGetArrayWrite(U, &u)); 159 d0 = uinit[0] - uinit[1]; 160 if (d0 == 0.0) q = ctx->k * t; 161 else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0; 162 u[0] = uinit[0] / (1.0 + uinit[1] * q); 163 u[1] = u[0] - d0; 164 u[2] = uinit[1] + uinit[2] - u[1]; 165 PetscCall(VecRestoreArrayWrite(U, &u)); 166 PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit)); 167 PetscFunctionReturn(PETSC_SUCCESS); 168 } 169 170 int main(int argc, char **argv) 171 { 172 TS ts; /* ODE integrator */ 173 Vec U; /* solution will be stored here */ 174 Mat A; /* Jacobian matrix */ 175 PetscMPIInt size; 176 PetscInt n = 3; 177 AppCtx ctx; 178 PetscScalar *u; 179 const char *const names[] = {"U1", "U2", "U3", NULL}; 180 PetscBool mf = PETSC_FALSE; 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 PetscCall(PetscOptionsGetBool(NULL, NULL, "-snes_mf_operator", &mf, NULL)); 191 192 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 193 Create necessary matrix and vectors 194 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 195 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 196 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE)); 197 PetscCall(MatSetFromOptions(A)); 198 PetscCall(MatSetUp(A)); 199 200 PetscCall(MatCreateVecs(A, &U, NULL)); 201 202 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 203 Set runtime options 204 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 205 ctx.k = .9; 206 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL)); 207 PetscCall(VecDuplicate(U, &ctx.initialsolution)); 208 PetscCall(VecGetArrayWrite(ctx.initialsolution, &u)); 209 u[0] = 1; 210 u[1] = .7; 211 u[2] = 0; 212 PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u)); 213 PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL)); 214 215 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 216 Create timestepping solver context 217 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 218 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 219 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 220 PetscCall(TSSetType(ts, TSROSW)); 221 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &ctx)); 222 if (!mf) PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx)); 223 else PetscCall(TSSetIJacobian(ts, NULL, NULL, (TSIJacobianFn *)IJacobian, &ctx)); 224 PetscCall(TSSetSolutionFunction(ts, (TSSolutionFn *)Solution, &ctx)); 225 226 { 227 DM dm; 228 void *ptr; 229 PetscCall(TSGetDM(ts, &dm)); 230 PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr)); 231 PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr)); 232 PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 233 PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 234 } 235 236 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 237 Set initial conditions 238 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 239 PetscCall(Solution(ts, 0, U, &ctx)); 240 PetscCall(TSSetSolution(ts, U)); 241 242 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243 Set solver options 244 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 245 PetscCall(TSSetTimeStep(ts, .001)); 246 PetscCall(TSSetMaxSteps(ts, 1000)); 247 PetscCall(TSSetMaxTime(ts, 20.0)); 248 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 249 PetscCall(TSSetFromOptions(ts)); 250 PetscCall(TSMonitorLGSetVariableNames(ts, names)); 251 252 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 253 Solve nonlinear system 254 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 255 PetscCall(TSSolve(ts, U)); 256 257 PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD)); 258 259 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 260 Free work space. All PETSc objects should be destroyed when they are no longer needed. 261 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 262 PetscCall(VecDestroy(&ctx.initialsolution)); 263 PetscCall(MatDestroy(&A)); 264 PetscCall(VecDestroy(&U)); 265 PetscCall(TSDestroy(&ts)); 266 267 PetscCall(PetscFinalize()); 268 return 0; 269 } 270 271 /*TEST 272 273 test: 274 args: -ts_view 275 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 276 277 test: 278 suffix: 2 279 args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 280 requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 281 output_file: output/ex1_1.out 282 283 test: 284 requires: !single 285 suffix: 3 286 args: -ts_view -snes_mf_operator 287 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 288 289 TEST*/ 290