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 /* 36c4762a1bSJed Brown Project: Generate a nicely formated 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 38*1baa6e33SBarry 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 78c4762a1bSJed Brown PetscErrorCode IFunctionView(AppCtx *ctx,PetscViewer v) 79c4762a1bSJed Brown { 80c4762a1bSJed Brown PetscFunctionBegin; 819566063dSJacob Faibussowitsch PetscCall(PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR)); 82c4762a1bSJed Brown PetscFunctionReturn(0); 83c4762a1bSJed Brown } 84c4762a1bSJed Brown 85c4762a1bSJed Brown PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v) 86c4762a1bSJed Brown { 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 */ 96c4762a1bSJed Brown PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 97c4762a1bSJed Brown { 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 */ 118c4762a1bSJed Brown PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx) 119c4762a1bSJed Brown { 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)); 127c4762a1bSJed Brown J[0][0] = a + ctx->k*u[1]; J[0][1] = ctx->k*u[0]; J[0][2] = 0.0; 128c4762a1bSJed Brown J[1][0] = ctx->k*u[1]; J[1][1] = a + ctx->k*u[0]; J[1][2] = 0.0; 129c4762a1bSJed Brown J[2][0] = -ctx->k*u[1]; J[2][1] = -ctx->k*u[0]; J[2][2] = a; 1309566063dSJacob Faibussowitsch PetscCall(MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES)); 1319566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(U,&u)); 1329566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(Udot,&udot)); 133c4762a1bSJed Brown 1349566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1359566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 136c4762a1bSJed Brown if (A != B) { 1379566063dSJacob Faibussowitsch PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1389566063dSJacob Faibussowitsch PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 139c4762a1bSJed Brown } 140c4762a1bSJed Brown PetscFunctionReturn(0); 141c4762a1bSJed Brown } 142c4762a1bSJed Brown 143c4762a1bSJed Brown /* 144c4762a1bSJed Brown Defines the exact (analytic) solution to the ODE 145c4762a1bSJed Brown */ 146c4762a1bSJed Brown static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx) 147c4762a1bSJed Brown { 148c4762a1bSJed Brown const PetscScalar *uinit; 149c4762a1bSJed Brown PetscScalar *u,d0,q; 150c4762a1bSJed Brown 151c4762a1bSJed Brown PetscFunctionBegin; 1529566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(ctx->initialsolution,&uinit)); 1539566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(U,&u)); 154c4762a1bSJed Brown d0 = uinit[0] - uinit[1]; 155c4762a1bSJed Brown if (d0 == 0.0) q = ctx->k*t; 156c4762a1bSJed Brown else q = (1.0 - PetscExpScalar(-ctx->k*t*d0))/d0; 157c4762a1bSJed Brown u[0] = uinit[0]/(1.0 + uinit[1]*q); 158c4762a1bSJed Brown u[1] = u[0] - d0; 159c4762a1bSJed Brown u[2] = uinit[1] + uinit[2] - u[1]; 1609566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(U,&u)); 1619566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(ctx->initialsolution,&uinit)); 162c4762a1bSJed Brown PetscFunctionReturn(0); 163c4762a1bSJed Brown } 164c4762a1bSJed Brown 165c4762a1bSJed Brown int main(int argc,char **argv) 166c4762a1bSJed Brown { 167c4762a1bSJed Brown TS ts; /* ODE integrator */ 168c4762a1bSJed Brown Vec U; /* solution will be stored here */ 169c4762a1bSJed Brown Mat A; /* Jacobian matrix */ 170c4762a1bSJed Brown PetscMPIInt size; 171c4762a1bSJed Brown PetscInt n = 3; 172c4762a1bSJed Brown AppCtx ctx; 173c4762a1bSJed Brown PetscScalar *u; 174c4762a1bSJed Brown const char * const names[] = {"U1","U2","U3",NULL}; 175c4762a1bSJed Brown 176c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 177c4762a1bSJed Brown Initialize program 178c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1799566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 1809566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 1813c633725SBarry Smith PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 182c4762a1bSJed Brown 183c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 184c4762a1bSJed Brown Create necessary matrix and vectors 185c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1869566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 1879566063dSJacob Faibussowitsch PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 1889566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(A)); 1899566063dSJacob Faibussowitsch PetscCall(MatSetUp(A)); 190c4762a1bSJed Brown 1919566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&U,NULL)); 192c4762a1bSJed Brown 193c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194c4762a1bSJed Brown Set runtime options 195c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 196c4762a1bSJed Brown ctx.k = .9; 1979566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 1989566063dSJacob Faibussowitsch PetscCall(VecDuplicate(U,&ctx.initialsolution)); 1999566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(ctx.initialsolution,&u)); 200c4762a1bSJed Brown u[0] = 1; 201c4762a1bSJed Brown u[1] = .7; 202c4762a1bSJed Brown u[2] = 0; 2039566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(ctx.initialsolution,&u)); 2049566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 205c4762a1bSJed Brown 206c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207c4762a1bSJed Brown Create timestepping solver context 208c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2099566063dSJacob Faibussowitsch PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 2109566063dSJacob Faibussowitsch PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 2119566063dSJacob Faibussowitsch PetscCall(TSSetType(ts,TSROSW)); 2129566063dSJacob Faibussowitsch PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 2139566063dSJacob Faibussowitsch PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 2149566063dSJacob Faibussowitsch PetscCall(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 215c4762a1bSJed Brown 216c4762a1bSJed Brown { 217c4762a1bSJed Brown DM dm; 218c4762a1bSJed Brown void *ptr; 2199566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts,&dm)); 2209566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL,"IFunctionView",&ptr)); 2219566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 2229566063dSJacob Faibussowitsch PetscCall(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 2239566063dSJacob Faibussowitsch PetscCall(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 224c4762a1bSJed Brown } 225c4762a1bSJed Brown 226c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 227c4762a1bSJed Brown Set initial conditions 228c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2299566063dSJacob Faibussowitsch PetscCall(Solution(ts,0,U,&ctx)); 2309566063dSJacob Faibussowitsch PetscCall(TSSetSolution(ts,U)); 231c4762a1bSJed Brown 232c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233c4762a1bSJed Brown Set solver options 234c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2359566063dSJacob Faibussowitsch PetscCall(TSSetTimeStep(ts,.001)); 2369566063dSJacob Faibussowitsch PetscCall(TSSetMaxSteps(ts,1000)); 2379566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts,20.0)); 2389566063dSJacob Faibussowitsch PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 2399566063dSJacob Faibussowitsch PetscCall(TSSetFromOptions(ts)); 2409566063dSJacob Faibussowitsch PetscCall(TSMonitorLGSetVariableNames(ts,names)); 241c4762a1bSJed Brown 242c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243c4762a1bSJed Brown Solve nonlinear system 244c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2459566063dSJacob Faibussowitsch PetscCall(TSSolve(ts,U)); 246c4762a1bSJed Brown 2479566063dSJacob Faibussowitsch PetscCall(TSView(ts,PETSC_VIEWER_BINARY_WORLD)); 248c4762a1bSJed Brown 249c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 250c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they are no longer needed. 251c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2529566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ctx.initialsolution)); 2539566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 2549566063dSJacob Faibussowitsch PetscCall(VecDestroy(&U)); 2559566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts)); 256c4762a1bSJed Brown 2579566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 258b122ec5aSJacob Faibussowitsch return 0; 259c4762a1bSJed Brown } 260c4762a1bSJed Brown 261c4762a1bSJed Brown /*TEST 262c4762a1bSJed Brown 263c4762a1bSJed Brown test: 264c4762a1bSJed Brown args: -ts_view 265dfd57a17SPierre Jolivet requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 266c4762a1bSJed Brown 267c4762a1bSJed Brown test: 268c4762a1bSJed Brown suffix: 2 269c4762a1bSJed Brown args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 270dfd57a17SPierre Jolivet requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 271c4762a1bSJed Brown output_file: output/ex1_1.out 272c4762a1bSJed Brown 273c4762a1bSJed Brown TEST*/ 274