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 PetscErrorCode ierr; 81c4762a1bSJed Brown 82c4762a1bSJed Brown PetscFunctionBegin; 83c4762a1bSJed Brown ierr = PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR);CHKERRQ(ierr); 84c4762a1bSJed Brown PetscFunctionReturn(0); 85c4762a1bSJed Brown } 86c4762a1bSJed Brown 87c4762a1bSJed Brown PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v) 88c4762a1bSJed Brown { 89c4762a1bSJed Brown PetscErrorCode ierr; 90c4762a1bSJed Brown 91c4762a1bSJed Brown PetscFunctionBegin; 92c4762a1bSJed Brown ierr = PetscNew(ctx);CHKERRQ(ierr); 93c4762a1bSJed Brown ierr = PetscViewerBinaryRead(v,&(*ctx)->k,1,NULL,PETSC_SCALAR);CHKERRQ(ierr); 94c4762a1bSJed Brown PetscFunctionReturn(0); 95c4762a1bSJed Brown } 96c4762a1bSJed Brown 97c4762a1bSJed Brown /* 98c4762a1bSJed Brown Defines the ODE passed to the ODE solver 99c4762a1bSJed Brown */ 100c4762a1bSJed Brown PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 101c4762a1bSJed Brown { 102c4762a1bSJed Brown PetscErrorCode ierr; 103c4762a1bSJed Brown PetscScalar *f; 104c4762a1bSJed Brown const PetscScalar *u,*udot; 105c4762a1bSJed Brown 106c4762a1bSJed Brown PetscFunctionBegin; 107c4762a1bSJed Brown /* The next three lines allow us to access the entries of the vectors directly */ 108c4762a1bSJed Brown ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 109c4762a1bSJed Brown ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 110303a5415SBarry Smith ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 111c4762a1bSJed Brown f[0] = udot[0] + ctx->k*u[0]*u[1]; 112c4762a1bSJed Brown f[1] = udot[1] + ctx->k*u[0]*u[1]; 113c4762a1bSJed Brown f[2] = udot[2] - ctx->k*u[0]*u[1]; 114c4762a1bSJed Brown ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 115c4762a1bSJed Brown ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 116303a5415SBarry Smith ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 117c4762a1bSJed Brown PetscFunctionReturn(0); 118c4762a1bSJed Brown } 119c4762a1bSJed Brown 120c4762a1bSJed Brown /* 121c4762a1bSJed Brown Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 122c4762a1bSJed Brown */ 123c4762a1bSJed Brown PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx) 124c4762a1bSJed Brown { 125c4762a1bSJed Brown PetscErrorCode ierr; 126c4762a1bSJed Brown PetscInt rowcol[] = {0,1,2}; 127c4762a1bSJed Brown PetscScalar J[3][3]; 128c4762a1bSJed Brown const PetscScalar *u,*udot; 129c4762a1bSJed Brown 130c4762a1bSJed Brown PetscFunctionBegin; 131c4762a1bSJed Brown ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 132c4762a1bSJed Brown ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 133c4762a1bSJed Brown J[0][0] = a + ctx->k*u[1]; J[0][1] = ctx->k*u[0]; J[0][2] = 0.0; 134c4762a1bSJed Brown J[1][0] = ctx->k*u[1]; J[1][1] = a + ctx->k*u[0]; J[1][2] = 0.0; 135c4762a1bSJed Brown J[2][0] = -ctx->k*u[1]; J[2][1] = -ctx->k*u[0]; J[2][2] = a; 136c4762a1bSJed Brown ierr = MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 137c4762a1bSJed Brown ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 138c4762a1bSJed Brown ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 139c4762a1bSJed Brown 140c4762a1bSJed Brown ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 141c4762a1bSJed Brown ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 142c4762a1bSJed Brown if (A != B) { 143c4762a1bSJed Brown ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 144c4762a1bSJed Brown ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 145c4762a1bSJed Brown } 146c4762a1bSJed Brown PetscFunctionReturn(0); 147c4762a1bSJed Brown } 148c4762a1bSJed Brown 149c4762a1bSJed Brown /* 150c4762a1bSJed Brown Defines the exact (analytic) solution to the ODE 151c4762a1bSJed Brown */ 152c4762a1bSJed Brown static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx) 153c4762a1bSJed Brown { 154c4762a1bSJed Brown PetscErrorCode ierr; 155c4762a1bSJed Brown const PetscScalar *uinit; 156c4762a1bSJed Brown PetscScalar *u,d0,q; 157c4762a1bSJed Brown 158c4762a1bSJed Brown PetscFunctionBegin; 159c4762a1bSJed Brown ierr = VecGetArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 160303a5415SBarry Smith ierr = VecGetArrayWrite(U,&u);CHKERRQ(ierr); 161c4762a1bSJed Brown d0 = uinit[0] - uinit[1]; 162c4762a1bSJed Brown if (d0 == 0.0) q = ctx->k*t; 163c4762a1bSJed Brown else q = (1.0 - PetscExpScalar(-ctx->k*t*d0))/d0; 164c4762a1bSJed Brown u[0] = uinit[0]/(1.0 + uinit[1]*q); 165c4762a1bSJed Brown u[1] = u[0] - d0; 166c4762a1bSJed Brown u[2] = uinit[1] + uinit[2] - u[1]; 167303a5415SBarry Smith ierr = VecRestoreArrayWrite(U,&u);CHKERRQ(ierr); 168c4762a1bSJed Brown ierr = VecRestoreArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 169c4762a1bSJed Brown PetscFunctionReturn(0); 170c4762a1bSJed Brown } 171c4762a1bSJed Brown 172c4762a1bSJed Brown int main(int argc,char **argv) 173c4762a1bSJed Brown { 174c4762a1bSJed Brown TS ts; /* ODE integrator */ 175c4762a1bSJed Brown Vec U; /* solution will be stored here */ 176c4762a1bSJed Brown Mat A; /* Jacobian matrix */ 177c4762a1bSJed Brown PetscErrorCode ierr; 178c4762a1bSJed Brown PetscMPIInt size; 179c4762a1bSJed Brown PetscInt n = 3; 180c4762a1bSJed Brown AppCtx ctx; 181c4762a1bSJed Brown PetscScalar *u; 182c4762a1bSJed Brown const char * const names[] = {"U1","U2","U3",NULL}; 183c4762a1bSJed Brown 184c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185c4762a1bSJed Brown Initialize program 186c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 188ffc4695bSBarry Smith ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 189c4762a1bSJed Brown if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 190c4762a1bSJed Brown 191c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 192c4762a1bSJed Brown Create necessary matrix and vectors 193c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 194c4762a1bSJed Brown ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 195c4762a1bSJed Brown ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 196c4762a1bSJed Brown ierr = MatSetFromOptions(A);CHKERRQ(ierr); 197c4762a1bSJed Brown ierr = MatSetUp(A);CHKERRQ(ierr); 198c4762a1bSJed Brown 199c4762a1bSJed Brown ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 200c4762a1bSJed Brown 201c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202c4762a1bSJed Brown Set runtime options 203c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 204c4762a1bSJed Brown ctx.k = .9; 205c4762a1bSJed Brown ierr = PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL);CHKERRQ(ierr); 206c4762a1bSJed Brown ierr = VecDuplicate(U,&ctx.initialsolution);CHKERRQ(ierr); 207303a5415SBarry Smith ierr = VecGetArrayWrite(ctx.initialsolution,&u);CHKERRQ(ierr); 208c4762a1bSJed Brown u[0] = 1; 209c4762a1bSJed Brown u[1] = .7; 210c4762a1bSJed Brown u[2] = 0; 211303a5415SBarry Smith ierr = VecRestoreArrayWrite(ctx.initialsolution,&u);CHKERRQ(ierr); 212c4762a1bSJed Brown ierr = PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL);CHKERRQ(ierr); 213c4762a1bSJed Brown 214c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215c4762a1bSJed Brown Create timestepping solver context 216c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217c4762a1bSJed Brown ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 218c4762a1bSJed Brown ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 219c4762a1bSJed Brown ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); 220c4762a1bSJed Brown ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx);CHKERRQ(ierr); 221c4762a1bSJed Brown ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);CHKERRQ(ierr); 222c4762a1bSJed Brown ierr = TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx);CHKERRQ(ierr); 223c4762a1bSJed Brown 224c4762a1bSJed Brown { 225c4762a1bSJed Brown DM dm; 226c4762a1bSJed Brown void *ptr; 227c4762a1bSJed Brown ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 228c4762a1bSJed Brown ierr = PetscDLSym(NULL,"IFunctionView",&ptr);CHKERRQ(ierr); 229c4762a1bSJed Brown ierr = PetscDLSym(NULL,"IFunctionLoad",&ptr);CHKERRQ(ierr); 230c4762a1bSJed Brown ierr = DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 231c4762a1bSJed Brown ierr = DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 232c4762a1bSJed Brown } 233c4762a1bSJed Brown 234c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 235c4762a1bSJed Brown Set initial conditions 236c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 237c4762a1bSJed Brown ierr = Solution(ts,0,U,&ctx);CHKERRQ(ierr); 238c4762a1bSJed Brown ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 239c4762a1bSJed Brown 240c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 241c4762a1bSJed Brown Set solver options 242c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 243c4762a1bSJed Brown ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr); 244c4762a1bSJed Brown ierr = TSSetMaxSteps(ts,1000);CHKERRQ(ierr); 245c4762a1bSJed Brown ierr = TSSetMaxTime(ts,20.0);CHKERRQ(ierr); 246c4762a1bSJed Brown ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 247c4762a1bSJed Brown ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 248c4762a1bSJed Brown ierr = TSMonitorLGSetVariableNames(ts,names);CHKERRQ(ierr); 249c4762a1bSJed Brown 250c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 251c4762a1bSJed Brown Solve nonlinear system 252c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 253c4762a1bSJed Brown ierr = TSSolve(ts,U);CHKERRQ(ierr); 254c4762a1bSJed Brown 255c4762a1bSJed Brown ierr = TSView(ts,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 256c4762a1bSJed Brown 257c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 258c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they are no longer needed. 259c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 260c4762a1bSJed Brown ierr = VecDestroy(&ctx.initialsolution);CHKERRQ(ierr); 261c4762a1bSJed Brown ierr = MatDestroy(&A);CHKERRQ(ierr); 262c4762a1bSJed Brown ierr = VecDestroy(&U);CHKERRQ(ierr); 263c4762a1bSJed Brown ierr = TSDestroy(&ts);CHKERRQ(ierr); 264c4762a1bSJed Brown 265c4762a1bSJed Brown ierr = PetscFinalize(); 266c4762a1bSJed Brown return ierr; 267c4762a1bSJed Brown } 268c4762a1bSJed Brown 269c4762a1bSJed Brown /*TEST 270c4762a1bSJed Brown 271c4762a1bSJed Brown test: 272c4762a1bSJed Brown args: -ts_view 273*dfd57a17SPierre Jolivet requires: dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 274c4762a1bSJed Brown 275c4762a1bSJed Brown test: 276c4762a1bSJed Brown suffix: 2 277c4762a1bSJed Brown args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 278*dfd57a17SPierre Jolivet requires: x dlsym defined(PETSC_HAVE_DYNAMIC_LIBRARIES) 279c4762a1bSJed Brown output_file: output/ex1_1.out 280c4762a1bSJed Brown 281c4762a1bSJed Brown TEST*/ 282