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 181 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 182 Initialize program 183 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 184 PetscFunctionBeginUser; 185 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 186 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 187 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 188 189 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 190 Create necessary matrix and vectors 191 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 192 PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 193 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE)); 194 PetscCall(MatSetFromOptions(A)); 195 PetscCall(MatSetUp(A)); 196 197 PetscCall(MatCreateVecs(A, &U, NULL)); 198 199 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 200 Set runtime options 201 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 202 ctx.k = .9; 203 PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL)); 204 PetscCall(VecDuplicate(U, &ctx.initialsolution)); 205 PetscCall(VecGetArrayWrite(ctx.initialsolution, &u)); 206 u[0] = 1; 207 u[1] = .7; 208 u[2] = 0; 209 PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u)); 210 PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL)); 211 212 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 213 Create timestepping solver context 214 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 215 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 216 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 217 PetscCall(TSSetType(ts, TSROSW)); 218 PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &ctx)); 219 PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobian)IJacobian, &ctx)); 220 PetscCall(TSSetSolutionFunction(ts, (TSSolutionFunction)Solution, &ctx)); 221 222 { 223 DM dm; 224 void *ptr; 225 PetscCall(TSGetDM(ts, &dm)); 226 PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr)); 227 PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr)); 228 PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 229 PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 230 } 231 232 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233 Set initial conditions 234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 235 PetscCall(Solution(ts, 0, U, &ctx)); 236 PetscCall(TSSetSolution(ts, U)); 237 238 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 239 Set solver options 240 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 241 PetscCall(TSSetTimeStep(ts, .001)); 242 PetscCall(TSSetMaxSteps(ts, 1000)); 243 PetscCall(TSSetMaxTime(ts, 20.0)); 244 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 245 PetscCall(TSSetFromOptions(ts)); 246 PetscCall(TSMonitorLGSetVariableNames(ts, names)); 247 248 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 249 Solve nonlinear system 250 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 251 PetscCall(TSSolve(ts, U)); 252 253 PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD)); 254 255 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 256 Free work space. All PETSc objects should be destroyed when they are no longer needed. 257 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 258 PetscCall(VecDestroy(&ctx.initialsolution)); 259 PetscCall(MatDestroy(&A)); 260 PetscCall(VecDestroy(&U)); 261 PetscCall(TSDestroy(&ts)); 262 263 PetscCall(PetscFinalize()); 264 return 0; 265 } 266 267 /*TEST 268 269 test: 270 args: -ts_view 271 requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 272 273 test: 274 suffix: 2 275 args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 276 requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 277 output_file: output/ex1_1.out 278 279 TEST*/ 280