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; 815f80ce2aSJacob 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; 885f80ce2aSJacob Faibussowitsch CHKERRQ(PetscNew(ctx)); 895f80ce2aSJacob 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 */ 1035f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(U,&u)); 1045f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(Udot,&udot)); 1055f80ce2aSJacob 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]; 1095f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(U,&u)); 1105f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 1115f80ce2aSJacob 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; 1255f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(U,&u)); 1265f80ce2aSJacob 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; 1305f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES)); 1315f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(U,&u)); 1325f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayRead(Udot,&udot)); 133c4762a1bSJed Brown 1345f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 1355f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 136c4762a1bSJed Brown if (A != B) { 1375f80ce2aSJacob Faibussowitsch CHKERRQ(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 1385f80ce2aSJacob 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; 1525f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayRead(ctx->initialsolution,&uinit)); 1535f80ce2aSJacob 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]; 1605f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayWrite(U,&u)); 1615f80ce2aSJacob 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 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*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscInitialize(&argc,&argv,(char*)0,help)); 1805f80ce2aSJacob Faibussowitsch CHKERRMPI(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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 1865f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreate(PETSC_COMM_WORLD,&A)); 1875f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 1885f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetFromOptions(A)); 1895f80ce2aSJacob Faibussowitsch CHKERRQ(MatSetUp(A)); 190c4762a1bSJed Brown 1915f80ce2aSJacob Faibussowitsch CHKERRQ(MatCreateVecs(A,&U,NULL)); 192c4762a1bSJed Brown 193c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 194c4762a1bSJed Brown Set runtime options 195c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 196c4762a1bSJed Brown ctx.k = .9; 1975f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL)); 1985f80ce2aSJacob Faibussowitsch CHKERRQ(VecDuplicate(U,&ctx.initialsolution)); 1995f80ce2aSJacob Faibussowitsch CHKERRQ(VecGetArrayWrite(ctx.initialsolution,&u)); 200c4762a1bSJed Brown u[0] = 1; 201c4762a1bSJed Brown u[1] = .7; 202c4762a1bSJed Brown u[2] = 0; 2035f80ce2aSJacob Faibussowitsch CHKERRQ(VecRestoreArrayWrite(ctx.initialsolution,&u)); 2045f80ce2aSJacob Faibussowitsch CHKERRQ(PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL)); 205c4762a1bSJed Brown 206c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 207c4762a1bSJed Brown Create timestepping solver context 208c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2095f80ce2aSJacob Faibussowitsch CHKERRQ(TSCreate(PETSC_COMM_WORLD,&ts)); 2105f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetProblemType(ts,TS_NONLINEAR)); 2115f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetType(ts,TSROSW)); 2125f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 2135f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 2145f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx)); 215c4762a1bSJed Brown 216c4762a1bSJed Brown { 217c4762a1bSJed Brown DM dm; 218c4762a1bSJed Brown void *ptr; 2195f80ce2aSJacob Faibussowitsch CHKERRQ(TSGetDM(ts,&dm)); 2205f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDLSym(NULL,"IFunctionView",&ptr)); 2215f80ce2aSJacob Faibussowitsch CHKERRQ(PetscDLSym(NULL,"IFunctionLoad",&ptr)); 2225f80ce2aSJacob Faibussowitsch CHKERRQ(DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 2235f80ce2aSJacob Faibussowitsch CHKERRQ(DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad)); 224c4762a1bSJed Brown } 225c4762a1bSJed Brown 226c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 227c4762a1bSJed Brown Set initial conditions 228c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2295f80ce2aSJacob Faibussowitsch CHKERRQ(Solution(ts,0,U,&ctx)); 2305f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetSolution(ts,U)); 231c4762a1bSJed Brown 232c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 233c4762a1bSJed Brown Set solver options 234c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2355f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetTimeStep(ts,.001)); 2365f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetMaxSteps(ts,1000)); 2375f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetMaxTime(ts,20.0)); 2385f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 2395f80ce2aSJacob Faibussowitsch CHKERRQ(TSSetFromOptions(ts)); 2405f80ce2aSJacob Faibussowitsch CHKERRQ(TSMonitorLGSetVariableNames(ts,names)); 241c4762a1bSJed Brown 242c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 243c4762a1bSJed Brown Solve nonlinear system 244c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2455f80ce2aSJacob Faibussowitsch CHKERRQ(TSSolve(ts,U)); 246c4762a1bSJed Brown 2475f80ce2aSJacob Faibussowitsch CHKERRQ(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 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2525f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&ctx.initialsolution)); 2535f80ce2aSJacob Faibussowitsch CHKERRQ(MatDestroy(&A)); 2545f80ce2aSJacob Faibussowitsch CHKERRQ(VecDestroy(&U)); 2555f80ce2aSJacob Faibussowitsch CHKERRQ(TSDestroy(&ts)); 256c4762a1bSJed Brown 257*b122ec5aSJacob Faibussowitsch CHKERRQ(PetscFinalize()); 258*b122ec5aSJacob 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