1c4762a1bSJed Brown static char help[] = "Nonlinear Reaction Problem from Chemistry.\n"; 2c4762a1bSJed Brown 3c4762a1bSJed Brown /*F 4c4762a1bSJed Brown 5c4762a1bSJed Brown This directory contains examples based on the PDES/ODES given in the book 6c4762a1bSJed Brown Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by 7c4762a1bSJed Brown W. Hundsdorf and J.G. Verwer 8c4762a1bSJed Brown 9c4762a1bSJed Brown Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry 10c4762a1bSJed Brown 11c4762a1bSJed Brown \begin{eqnarray} 12c4762a1bSJed Brown {U_1}_t - k U_1 U_2 & = & 0 \\ 13c4762a1bSJed Brown {U_2}_t - k U_1 U_2 & = & 0 \\ 14c4762a1bSJed Brown {U_3}_t - k U_1 U_2 & = & 0 15c4762a1bSJed Brown \end{eqnarray} 16c4762a1bSJed Brown 17c4762a1bSJed Brown Helpful runtime monitoring options: 18c4762a1bSJed Brown -ts_view - prints information about the solver being used 19da81f932SPierre Jolivet -ts_monitor - prints the progress of the solver 20c4762a1bSJed Brown -ts_adapt_monitor - prints the progress of the time-step adaptor 21c4762a1bSJed Brown -ts_monitor_lg_timestep - plots the size of each timestep (at each time-step) 22c4762a1bSJed Brown -ts_monitor_lg_solution - plots each component of the solution as a function of time (at each timestep) 23c4762a1bSJed Brown -ts_monitor_lg_error - plots each component of the error in the solution as a function of time (at each timestep) 24c4762a1bSJed Brown -draw_pause -2 - hold the plots a the end of the solution process, enter a mouse press in each window to end the process 25c4762a1bSJed Brown 26c4762a1bSJed Brown -ts_monitor_lg_timestep -1 - plots the size of each timestep (at the end of the solution process) 27c4762a1bSJed Brown -ts_monitor_lg_solution -1 - plots each component of the solution as a function of time (at the end of the solution process) 28c4762a1bSJed Brown -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) 29c4762a1bSJed Brown -lg_use_markers false - do NOT show the data points on the plots 30c4762a1bSJed Brown -draw_save - save the timestep and solution plot as a .Gif image file 31c4762a1bSJed Brown 32c4762a1bSJed Brown F*/ 33c4762a1bSJed Brown 34c4762a1bSJed Brown /* 3535cb6cd3SPierre Jolivet Project: Generate a nicely formatted HTML page using 36c4762a1bSJed Brown 1) the HTML version of the source code and text in this file, use make html to generate the file ex1.c.html 371baa6e33SBarry Smith 2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_$_1_0.Gif) and 38c4762a1bSJed Brown 3) the text output (output.txt) generated by running the following commands. 39c4762a1bSJed Brown 4) <iframe src="generated_topics.html" scrolling="no" frameborder="0" width=600 height=300></iframe> 40c4762a1bSJed Brown 41c4762a1bSJed Brown rm -rf *.Gif 42c4762a1bSJed Brown ./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 43c4762a1bSJed Brown 44c4762a1bSJed Brown For example something like 45c4762a1bSJed Brown <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> 46c4762a1bSJed Brown <html> 47c4762a1bSJed Brown <head> 48c4762a1bSJed Brown <meta http-equiv="content-type" content="text/html;charset=utf-8"> 49c4762a1bSJed Brown <title>PETSc Example -- Nonlinear Reaction Problem from Chemistry</title> 50c4762a1bSJed Brown </head> 51c4762a1bSJed Brown <body> 52c4762a1bSJed Brown <iframe src="ex1.c.html" scrolling="yes" frameborder="1" width=2000 height=400></iframe> 53c4762a1bSJed Brown <img alt="" src="Draw_0x84000000_0_0.Gif"/><img alt="" src="Draw_0x84000001_0_0.Gif"/><img alt="" src="Draw_0x84000001_1_0.Gif"/> 54c4762a1bSJed Brown <iframe src="output.txt" scrolling="yes" frameborder="1" width=2000 height=1000></iframe> 55c4762a1bSJed Brown </body> 56c4762a1bSJed Brown </html> 57c4762a1bSJed Brown 58c4762a1bSJed Brown */ 59c4762a1bSJed Brown 60c4762a1bSJed Brown /* 61c4762a1bSJed Brown Include "petscts.h" so that we can use TS solvers. Note that this 62c4762a1bSJed Brown file automatically includes: 63c4762a1bSJed Brown petscsys.h - base PETSc routines petscvec.h - vectors 64c4762a1bSJed Brown petscmat.h - matrices 65c4762a1bSJed Brown petscis.h - index sets petscksp.h - Krylov subspace methods 66c4762a1bSJed Brown petscviewer.h - viewers petscpc.h - preconditioners 67c4762a1bSJed Brown petscksp.h - linear solvers 68c4762a1bSJed Brown */ 69c4762a1bSJed Brown 70c4762a1bSJed Brown #include <petscts.h> 71c4762a1bSJed Brown 72c4762a1bSJed Brown typedef struct { 73c4762a1bSJed Brown PetscScalar k; 74c4762a1bSJed Brown Vec initialsolution; 75c4762a1bSJed Brown } AppCtx; 76c4762a1bSJed Brown 77d71ae5a4SJacob Faibussowitsch PetscErrorCode IFunctionView(AppCtx *ctx, PetscViewer v) 78d71ae5a4SJacob Faibussowitsch { 79c4762a1bSJed Brown PetscFunctionBegin; 809566063dSJacob Faibussowitsch PetscCall(PetscViewerBinaryWrite(v, &ctx->k, 1, PETSC_SCALAR)); 813ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 82c4762a1bSJed Brown } 83c4762a1bSJed Brown 84d71ae5a4SJacob Faibussowitsch PetscErrorCode IFunctionLoad(AppCtx **ctx, PetscViewer v) 85d71ae5a4SJacob Faibussowitsch { 86c4762a1bSJed Brown PetscFunctionBegin; 879566063dSJacob Faibussowitsch PetscCall(PetscNew(ctx)); 889566063dSJacob Faibussowitsch PetscCall(PetscViewerBinaryRead(v, &(*ctx)->k, 1, NULL, PETSC_SCALAR)); 893ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 90c4762a1bSJed Brown } 91c4762a1bSJed Brown 92c4762a1bSJed Brown /* 93c4762a1bSJed Brown Defines the ODE passed to the ODE solver 94c4762a1bSJed Brown */ 95d71ae5a4SJacob Faibussowitsch PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, AppCtx *ctx) 96d71ae5a4SJacob Faibussowitsch { 97c4762a1bSJed Brown PetscScalar *f; 98c4762a1bSJed Brown const PetscScalar *u, *udot; 99c4762a1bSJed Brown 100c4762a1bSJed Brown PetscFunctionBegin; 101c4762a1bSJed Brown /* The next three lines allow us to access the entries of the vectors directly */ 1029566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(U, &u)); 1039566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Udot, &udot)); 1049566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(F, &f)); 105c4762a1bSJed Brown f[0] = udot[0] + ctx->k * u[0] * u[1]; 106c4762a1bSJed Brown f[1] = udot[1] + ctx->k * u[0] * u[1]; 107c4762a1bSJed Brown f[2] = udot[2] - ctx->k * u[0] * u[1]; 1089566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(U, &u)); 1099566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Udot, &udot)); 1109566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(F, &f)); 1113ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 112c4762a1bSJed Brown } 113c4762a1bSJed Brown 114c4762a1bSJed Brown /* 115c4762a1bSJed Brown Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 116c4762a1bSJed Brown */ 117d71ae5a4SJacob Faibussowitsch PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx) 118d71ae5a4SJacob Faibussowitsch { 119c4762a1bSJed Brown PetscInt rowcol[] = {0, 1, 2}; 120c4762a1bSJed Brown PetscScalar J[3][3]; 121c4762a1bSJed Brown const PetscScalar *u, *udot; 122c4762a1bSJed Brown 123c4762a1bSJed Brown PetscFunctionBegin; 1249566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(U, &u)); 1259566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(Udot, &udot)); 1269371c9d4SSatish Balay J[0][0] = a + ctx->k * u[1]; 1279371c9d4SSatish Balay J[0][1] = ctx->k * u[0]; 1289371c9d4SSatish Balay J[0][2] = 0.0; 1299371c9d4SSatish Balay J[1][0] = ctx->k * u[1]; 1309371c9d4SSatish Balay J[1][1] = a + ctx->k * u[0]; 1319371c9d4SSatish Balay J[1][2] = 0.0; 1329371c9d4SSatish Balay J[2][0] = -ctx->k * u[1]; 1339371c9d4SSatish Balay J[2][1] = -ctx->k * u[0]; 1349371c9d4SSatish Balay J[2][2] = a; 1359566063dSJacob Faibussowitsch PetscCall(MatSetValues(B, 3, rowcol, 3, rowcol, &J[0][0], INSERT_VALUES)); 1369566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(U, &u)); 1379566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Udot, &udot)); 138c4762a1bSJed Brown 1399566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY)); 1409566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY)); 141c4762a1bSJed Brown if (A != B) { 1429566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY)); 1439566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY)); 144c4762a1bSJed Brown } 1453ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 146c4762a1bSJed Brown } 147c4762a1bSJed Brown 148c4762a1bSJed Brown /* 149c4762a1bSJed Brown Defines the exact (analytic) solution to the ODE 150c4762a1bSJed Brown */ 151d71ae5a4SJacob Faibussowitsch static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx) 152d71ae5a4SJacob Faibussowitsch { 153c4762a1bSJed Brown const PetscScalar *uinit; 154c4762a1bSJed Brown PetscScalar *u, d0, q; 155c4762a1bSJed Brown 156c4762a1bSJed Brown PetscFunctionBegin; 1579566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(ctx->initialsolution, &uinit)); 1589566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(U, &u)); 159c4762a1bSJed Brown d0 = uinit[0] - uinit[1]; 160c4762a1bSJed Brown if (d0 == 0.0) q = ctx->k * t; 161c4762a1bSJed Brown else q = (1.0 - PetscExpScalar(-ctx->k * t * d0)) / d0; 162c4762a1bSJed Brown u[0] = uinit[0] / (1.0 + uinit[1] * q); 163c4762a1bSJed Brown u[1] = u[0] - d0; 164c4762a1bSJed Brown u[2] = uinit[1] + uinit[2] - u[1]; 1659566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(U, &u)); 1669566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(ctx->initialsolution, &uinit)); 1673ba16761SJacob Faibussowitsch PetscFunctionReturn(PETSC_SUCCESS); 168c4762a1bSJed Brown } 169c4762a1bSJed Brown 170d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv) 171d71ae5a4SJacob Faibussowitsch { 172c4762a1bSJed Brown TS ts; /* ODE integrator */ 173c4762a1bSJed Brown Vec U; /* solution will be stored here */ 174c4762a1bSJed Brown Mat A; /* Jacobian matrix */ 175c4762a1bSJed Brown PetscMPIInt size; 176c4762a1bSJed Brown PetscInt n = 3; 177c4762a1bSJed Brown AppCtx ctx; 178c4762a1bSJed Brown PetscScalar *u; 179c4762a1bSJed Brown const char *const names[] = {"U1", "U2", "U3", NULL}; 180c4762a1bSJed Brown 181c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 182c4762a1bSJed Brown Initialize program 183c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 184327415f7SBarry Smith PetscFunctionBeginUser; 1859566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 1869566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size)); 1873c633725SBarry Smith PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs"); 188c4762a1bSJed Brown 189c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 190c4762a1bSJed Brown Create necessary matrix and vectors 191c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1929566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD, &A)); 1939566063dSJacob Faibussowitsch PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE)); 1949566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(A)); 1959566063dSJacob Faibussowitsch PetscCall(MatSetUp(A)); 196c4762a1bSJed Brown 1979566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A, &U, NULL)); 198c4762a1bSJed Brown 199c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 200c4762a1bSJed Brown Set runtime options 201c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 202c4762a1bSJed Brown ctx.k = .9; 2039566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetScalar(NULL, NULL, "-k", &ctx.k, NULL)); 2049566063dSJacob Faibussowitsch PetscCall(VecDuplicate(U, &ctx.initialsolution)); 2059566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(ctx.initialsolution, &u)); 206c4762a1bSJed Brown u[0] = 1; 207c4762a1bSJed Brown u[1] = .7; 208c4762a1bSJed Brown u[2] = 0; 2099566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u)); 2109566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetVec(NULL, NULL, "-initial", ctx.initialsolution, NULL)); 211c4762a1bSJed Brown 212c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 213c4762a1bSJed Brown Create timestepping solver context 214c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2159566063dSJacob Faibussowitsch PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 2169566063dSJacob Faibussowitsch PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 2179566063dSJacob Faibussowitsch PetscCall(TSSetType(ts, TSROSW)); 218*8434afd1SBarry Smith PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &ctx)); 219*8434afd1SBarry Smith PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &ctx)); 220*8434afd1SBarry Smith PetscCall(TSSetSolutionFunction(ts, (TSSolutionFn *)Solution, &ctx)); 221c4762a1bSJed Brown 222c4762a1bSJed Brown { 223c4762a1bSJed Brown DM dm; 224c4762a1bSJed Brown void *ptr; 2259566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts, &dm)); 2269566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL, "IFunctionView", &ptr)); 2279566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL, "IFunctionLoad", &ptr)); 2289566063dSJacob Faibussowitsch PetscCall(DMTSSetIFunctionSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 2299566063dSJacob Faibussowitsch PetscCall(DMTSSetIJacobianSerialize(dm, (PetscErrorCode(*)(void *, PetscViewer))IFunctionView, (PetscErrorCode(*)(void **, PetscViewer))IFunctionLoad)); 230c4762a1bSJed Brown } 231c4762a1bSJed Brown 232c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233c4762a1bSJed Brown Set initial conditions 234c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2359566063dSJacob Faibussowitsch PetscCall(Solution(ts, 0, U, &ctx)); 2369566063dSJacob Faibussowitsch PetscCall(TSSetSolution(ts, U)); 237c4762a1bSJed Brown 238c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 239c4762a1bSJed Brown Set solver options 240c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2419566063dSJacob Faibussowitsch PetscCall(TSSetTimeStep(ts, .001)); 2429566063dSJacob Faibussowitsch PetscCall(TSSetMaxSteps(ts, 1000)); 2439566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts, 20.0)); 2449566063dSJacob Faibussowitsch PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 2459566063dSJacob Faibussowitsch PetscCall(TSSetFromOptions(ts)); 2469566063dSJacob Faibussowitsch PetscCall(TSMonitorLGSetVariableNames(ts, names)); 247c4762a1bSJed Brown 248c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 249c4762a1bSJed Brown Solve nonlinear system 250c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2519566063dSJacob Faibussowitsch PetscCall(TSSolve(ts, U)); 252c4762a1bSJed Brown 2539566063dSJacob Faibussowitsch PetscCall(TSView(ts, PETSC_VIEWER_BINARY_WORLD)); 254c4762a1bSJed Brown 255c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 256c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they are no longer needed. 257c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2589566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ctx.initialsolution)); 2599566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 2609566063dSJacob Faibussowitsch PetscCall(VecDestroy(&U)); 2619566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts)); 262c4762a1bSJed Brown 2639566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 264b122ec5aSJacob Faibussowitsch return 0; 265c4762a1bSJed Brown } 266c4762a1bSJed Brown 267c4762a1bSJed Brown /*TEST 268c4762a1bSJed Brown 269c4762a1bSJed Brown test: 270c4762a1bSJed Brown args: -ts_view 271dfd57a17SPierre Jolivet requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 272c4762a1bSJed Brown 273c4762a1bSJed Brown test: 274c4762a1bSJed Brown suffix: 2 275c4762a1bSJed Brown args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 276dfd57a17SPierre Jolivet requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 277c4762a1bSJed Brown output_file: output/ex1_1.out 278c4762a1bSJed Brown 279c4762a1bSJed Brown TEST*/ 280