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 38c4762a1bSJed Brown 2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_ZZZ_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; 81*5f80ce2aSJacob Faibussowitsch CHKERRQ(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; 88*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscNew(ctx)); 89*5f80ce2aSJacob Faibussowitsch CHKERRQ(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 */ 103*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(U,&u)); 104*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Udot,&udot)); 105*5f80ce2aSJacob Faibussowitsch CHKERRQ(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]; 109*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(U,&u)); 110*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 111*5f80ce2aSJacob Faibussowitsch CHKERRQ(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; 125*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(U,&u)); 126*5f80ce2aSJacob Faibussowitsch CHKERRQ(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; 130*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES)); 131*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(U,&u)); 132*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 133c4762a1bSJed Brown 134*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 135*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 136c4762a1bSJed Brown if (A != B) { 137*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 138*5f80ce2aSJacob Faibussowitsch CHKERRQ(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; 152*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(ctx->initialsolution,&uinit)); 153*5f80ce2aSJacob Faibussowitsch CHKERRQ(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]; 160*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayWrite(U,&u)); 161*5f80ce2aSJacob Faibussowitsch CHKERRQ(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 PetscErrorCode ierr; 171c4762a1bSJed Brown PetscMPIInt size; 172c4762a1bSJed Brown PetscInt n = 3; 173c4762a1bSJed Brown AppCtx ctx; 174c4762a1bSJed Brown PetscScalar *u; 175c4762a1bSJed Brown const char * const names[] = {"U1","U2","U3",NULL}; 176c4762a1bSJed Brown 177c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 178c4762a1bSJed Brown Initialize program 179c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 180c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 181*5f80ce2aSJacob Faibussowitsch CHKERRMPI(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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 188*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 189*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(A)); 190*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(A)); 191c4762a1bSJed Brown 192*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateVecs(A,&U,NULL)); 193c4762a1bSJed Brown 194c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 195c4762a1bSJed Brown Set runtime options 196c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 197c4762a1bSJed Brown ctx.k = .9; 198*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 199*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDuplicate(U,&ctx.initialsolution)); 200*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayWrite(ctx.initialsolution,&u)); 201c4762a1bSJed Brown u[0] = 1; 202c4762a1bSJed Brown u[1] = .7; 203c4762a1bSJed Brown u[2] = 0; 204*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayWrite(ctx.initialsolution,&u)); 205*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 206c4762a1bSJed Brown 207c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 208c4762a1bSJed Brown Create timestepping solver context 209c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 210*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 211*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 212*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetType(ts,TSROSW)); 213*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 214*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 215*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 216c4762a1bSJed Brown 217c4762a1bSJed Brown { 218c4762a1bSJed Brown DM dm; 219c4762a1bSJed Brown void *ptr; 220*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSGetDM(ts,&dm)); 221*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDLSym(NULL,"IFunctionView",&ptr)); 222*5f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 223*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 224*5f80ce2aSJacob Faibussowitsch CHKERRQ(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 225c4762a1bSJed Brown } 226c4762a1bSJed Brown 227c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 228c4762a1bSJed Brown Set initial conditions 229c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 230*5f80ce2aSJacob Faibussowitsch CHKERRQ(Solution(ts,0,U,&ctx)); 231*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetSolution(ts,U)); 232c4762a1bSJed Brown 233c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 234c4762a1bSJed Brown Set solver options 235c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 236*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetTimeStep(ts,.001)); 237*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetMaxSteps(ts,1000)); 238*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetMaxTime(ts,20.0)); 239*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 240*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetFromOptions(ts)); 241*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSMonitorLGSetVariableNames(ts,names)); 242c4762a1bSJed Brown 243c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 244c4762a1bSJed Brown Solve nonlinear system 245c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 246*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSSolve(ts,U)); 247c4762a1bSJed Brown 248*5f80ce2aSJacob Faibussowitsch CHKERRQ(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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 253*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&ctx.initialsolution)); 254*5f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A)); 255*5f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&U)); 256*5f80ce2aSJacob Faibussowitsch CHKERRQ(TSDestroy(&ts)); 257c4762a1bSJed Brown 258c4762a1bSJed Brown ierr = PetscFinalize(); 259c4762a1bSJed Brown return ierr; 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