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 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 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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 179*327415f7SBarry Smith PetscFunctionBeginUser; 1809566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 1819566063dSJacob Faibussowitsch PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 1823c633725SBarry Smith PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 183c4762a1bSJed Brown 184c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185c4762a1bSJed Brown Create necessary matrix and vectors 186c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1879566063dSJacob Faibussowitsch PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 1889566063dSJacob Faibussowitsch PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 1899566063dSJacob Faibussowitsch PetscCall(MatSetFromOptions(A)); 1909566063dSJacob Faibussowitsch PetscCall(MatSetUp(A)); 191c4762a1bSJed Brown 1929566063dSJacob Faibussowitsch PetscCall(MatCreateVecs(A,&U,NULL)); 193c4762a1bSJed Brown 194c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195c4762a1bSJed Brown Set runtime options 196c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197c4762a1bSJed Brown ctx.k = .9; 1989566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 1999566063dSJacob Faibussowitsch PetscCall(VecDuplicate(U,&ctx.initialsolution)); 2009566063dSJacob Faibussowitsch PetscCall(VecGetArrayWrite(ctx.initialsolution,&u)); 201c4762a1bSJed Brown u[0] = 1; 202c4762a1bSJed Brown u[1] = .7; 203c4762a1bSJed Brown u[2] = 0; 2049566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayWrite(ctx.initialsolution,&u)); 2059566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 206c4762a1bSJed Brown 207c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 208c4762a1bSJed Brown Create timestepping solver context 209c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2109566063dSJacob Faibussowitsch PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 2119566063dSJacob Faibussowitsch PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 2129566063dSJacob Faibussowitsch PetscCall(TSSetType(ts,TSROSW)); 2139566063dSJacob Faibussowitsch PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 2149566063dSJacob Faibussowitsch PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 2159566063dSJacob Faibussowitsch PetscCall(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 216c4762a1bSJed Brown 217c4762a1bSJed Brown { 218c4762a1bSJed Brown DM dm; 219c4762a1bSJed Brown void *ptr; 2209566063dSJacob Faibussowitsch PetscCall(TSGetDM(ts,&dm)); 2219566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL,"IFunctionView",&ptr)); 2229566063dSJacob Faibussowitsch PetscCall(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 2239566063dSJacob Faibussowitsch PetscCall(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 2249566063dSJacob Faibussowitsch PetscCall(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 225c4762a1bSJed Brown } 226c4762a1bSJed Brown 227c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 228c4762a1bSJed Brown Set initial conditions 229c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2309566063dSJacob Faibussowitsch PetscCall(Solution(ts,0,U,&ctx)); 2319566063dSJacob Faibussowitsch PetscCall(TSSetSolution(ts,U)); 232c4762a1bSJed Brown 233c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 234c4762a1bSJed Brown Set solver options 235c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2369566063dSJacob Faibussowitsch PetscCall(TSSetTimeStep(ts,.001)); 2379566063dSJacob Faibussowitsch PetscCall(TSSetMaxSteps(ts,1000)); 2389566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts,20.0)); 2399566063dSJacob Faibussowitsch PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 2409566063dSJacob Faibussowitsch PetscCall(TSSetFromOptions(ts)); 2419566063dSJacob Faibussowitsch PetscCall(TSMonitorLGSetVariableNames(ts,names)); 242c4762a1bSJed Brown 243c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244c4762a1bSJed Brown Solve nonlinear system 245c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2469566063dSJacob Faibussowitsch PetscCall(TSSolve(ts,U)); 247c4762a1bSJed Brown 2489566063dSJacob Faibussowitsch PetscCall(TSView(ts,PETSC_VIEWER_BINARY_WORLD)); 249c4762a1bSJed Brown 250c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 251c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they are no longer needed. 252c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2539566063dSJacob Faibussowitsch PetscCall(VecDestroy(&ctx.initialsolution)); 2549566063dSJacob Faibussowitsch PetscCall(MatDestroy(&A)); 2559566063dSJacob Faibussowitsch PetscCall(VecDestroy(&U)); 2569566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts)); 257c4762a1bSJed Brown 2589566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 259b122ec5aSJacob Faibussowitsch return 0; 260c4762a1bSJed Brown } 261c4762a1bSJed Brown 262c4762a1bSJed Brown /*TEST 263c4762a1bSJed Brown 264c4762a1bSJed Brown test: 265c4762a1bSJed Brown args: -ts_view 266dfd57a17SPierre Jolivet requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 267c4762a1bSJed Brown 268c4762a1bSJed Brown test: 269c4762a1bSJed Brown suffix: 2 270c4762a1bSJed Brown args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 271dfd57a17SPierre Jolivet requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 272c4762a1bSJed Brown output_file: output/ex1_1.out 273c4762a1bSJed Brown 274c4762a1bSJed Brown TEST*/ 275