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