1*c4762a1bSJed Brown 2*c4762a1bSJed Brown static char help[] = "Nonlinear Reaction Problem from Chemistry.\n"; 3*c4762a1bSJed Brown 4*c4762a1bSJed Brown /*F 5*c4762a1bSJed Brown 6*c4762a1bSJed Brown This directory contains examples based on the PDES/ODES given in the book 7*c4762a1bSJed Brown Numerical Solution of Time-Dependent Advection-Diffusion-Reaction Equations by 8*c4762a1bSJed Brown W. Hundsdorf and J.G. Verwer 9*c4762a1bSJed Brown 10*c4762a1bSJed Brown Page 3, Section 1.1 Nonlinear Reaction Problems from Chemistry 11*c4762a1bSJed Brown 12*c4762a1bSJed Brown \begin{eqnarray} 13*c4762a1bSJed Brown {U_1}_t - k U_1 U_2 & = & 0 \\ 14*c4762a1bSJed Brown {U_2}_t - k U_1 U_2 & = & 0 \\ 15*c4762a1bSJed Brown {U_3}_t - k U_1 U_2 & = & 0 16*c4762a1bSJed Brown \end{eqnarray} 17*c4762a1bSJed Brown 18*c4762a1bSJed Brown Helpful runtime monitoring options: 19*c4762a1bSJed Brown -ts_view - prints information about the solver being used 20*c4762a1bSJed Brown -ts_monitor - prints the progess of the solver 21*c4762a1bSJed Brown -ts_adapt_monitor - prints the progress of the time-step adaptor 22*c4762a1bSJed Brown -ts_monitor_lg_timestep - plots the size of each timestep (at each time-step) 23*c4762a1bSJed Brown -ts_monitor_lg_solution - plots each component of the solution as a function of time (at each timestep) 24*c4762a1bSJed Brown -ts_monitor_lg_error - plots each component of the error in the solution as a function of time (at each timestep) 25*c4762a1bSJed 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 26*c4762a1bSJed Brown 27*c4762a1bSJed Brown -ts_monitor_lg_timestep -1 - plots the size of each timestep (at the end of the solution process) 28*c4762a1bSJed Brown -ts_monitor_lg_solution -1 - plots each component of the solution as a function of time (at the end of the solution process) 29*c4762a1bSJed 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) 30*c4762a1bSJed Brown -lg_use_markers false - do NOT show the data points on the plots 31*c4762a1bSJed Brown -draw_save - save the timestep and solution plot as a .Gif image file 32*c4762a1bSJed Brown 33*c4762a1bSJed Brown F*/ 34*c4762a1bSJed Brown 35*c4762a1bSJed Brown /* 36*c4762a1bSJed Brown Project: Generate a nicely formated HTML page using 37*c4762a1bSJed 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*c4762a1bSJed Brown 2) the images (Draw_XXX_0_0.Gif Draw_YYY_0_0.Gif Draw_ZZZ_1_0.Gif) and 39*c4762a1bSJed Brown 3) the text output (output.txt) generated by running the following commands. 40*c4762a1bSJed Brown 4) <iframe src="generated_topics.html" scrolling="no" frameborder="0" width=600 height=300></iframe> 41*c4762a1bSJed Brown 42*c4762a1bSJed Brown rm -rf *.Gif 43*c4762a1bSJed 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 44*c4762a1bSJed Brown 45*c4762a1bSJed Brown For example something like 46*c4762a1bSJed Brown <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN"> 47*c4762a1bSJed Brown <html> 48*c4762a1bSJed Brown <head> 49*c4762a1bSJed Brown <meta http-equiv="content-type" content="text/html;charset=utf-8"> 50*c4762a1bSJed Brown <title>PETSc Example -- Nonlinear Reaction Problem from Chemistry</title> 51*c4762a1bSJed Brown </head> 52*c4762a1bSJed Brown <body> 53*c4762a1bSJed Brown <iframe src="ex1.c.html" scrolling="yes" frameborder="1" width=2000 height=400></iframe> 54*c4762a1bSJed 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"/> 55*c4762a1bSJed Brown <iframe src="output.txt" scrolling="yes" frameborder="1" width=2000 height=1000></iframe> 56*c4762a1bSJed Brown </body> 57*c4762a1bSJed Brown </html> 58*c4762a1bSJed Brown 59*c4762a1bSJed Brown */ 60*c4762a1bSJed Brown 61*c4762a1bSJed Brown /* 62*c4762a1bSJed Brown Include "petscts.h" so that we can use TS solvers. Note that this 63*c4762a1bSJed Brown file automatically includes: 64*c4762a1bSJed Brown petscsys.h - base PETSc routines petscvec.h - vectors 65*c4762a1bSJed Brown petscmat.h - matrices 66*c4762a1bSJed Brown petscis.h - index sets petscksp.h - Krylov subspace methods 67*c4762a1bSJed Brown petscviewer.h - viewers petscpc.h - preconditioners 68*c4762a1bSJed Brown petscksp.h - linear solvers 69*c4762a1bSJed Brown */ 70*c4762a1bSJed Brown 71*c4762a1bSJed Brown #include <petscts.h> 72*c4762a1bSJed Brown 73*c4762a1bSJed Brown typedef struct { 74*c4762a1bSJed Brown PetscScalar k; 75*c4762a1bSJed Brown Vec initialsolution; 76*c4762a1bSJed Brown } AppCtx; 77*c4762a1bSJed Brown 78*c4762a1bSJed Brown PetscErrorCode IFunctionView(AppCtx *ctx,PetscViewer v) 79*c4762a1bSJed Brown { 80*c4762a1bSJed Brown PetscErrorCode ierr; 81*c4762a1bSJed Brown 82*c4762a1bSJed Brown PetscFunctionBegin; 83*c4762a1bSJed Brown ierr = PetscViewerBinaryWrite(v,&ctx->k,1,PETSC_SCALAR);CHKERRQ(ierr); 84*c4762a1bSJed Brown PetscFunctionReturn(0); 85*c4762a1bSJed Brown } 86*c4762a1bSJed Brown 87*c4762a1bSJed Brown PetscErrorCode IFunctionLoad(AppCtx **ctx,PetscViewer v) 88*c4762a1bSJed Brown { 89*c4762a1bSJed Brown PetscErrorCode ierr; 90*c4762a1bSJed Brown 91*c4762a1bSJed Brown PetscFunctionBegin; 92*c4762a1bSJed Brown ierr = PetscNew(ctx);CHKERRQ(ierr); 93*c4762a1bSJed Brown ierr = PetscViewerBinaryRead(v,&(*ctx)->k,1,NULL,PETSC_SCALAR);CHKERRQ(ierr); 94*c4762a1bSJed Brown PetscFunctionReturn(0); 95*c4762a1bSJed Brown } 96*c4762a1bSJed Brown 97*c4762a1bSJed Brown /* 98*c4762a1bSJed Brown Defines the ODE passed to the ODE solver 99*c4762a1bSJed Brown */ 100*c4762a1bSJed Brown PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 101*c4762a1bSJed Brown { 102*c4762a1bSJed Brown PetscErrorCode ierr; 103*c4762a1bSJed Brown PetscScalar *f; 104*c4762a1bSJed Brown const PetscScalar *u,*udot; 105*c4762a1bSJed Brown 106*c4762a1bSJed Brown PetscFunctionBegin; 107*c4762a1bSJed Brown /* The next three lines allow us to access the entries of the vectors directly */ 108*c4762a1bSJed Brown ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 109*c4762a1bSJed Brown ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 110*c4762a1bSJed Brown ierr = VecGetArray(F,&f);CHKERRQ(ierr); 111*c4762a1bSJed Brown f[0] = udot[0] + ctx->k*u[0]*u[1]; 112*c4762a1bSJed Brown f[1] = udot[1] + ctx->k*u[0]*u[1]; 113*c4762a1bSJed Brown f[2] = udot[2] - ctx->k*u[0]*u[1]; 114*c4762a1bSJed Brown ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 115*c4762a1bSJed Brown ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 116*c4762a1bSJed Brown ierr = VecRestoreArray(F,&f);CHKERRQ(ierr); 117*c4762a1bSJed Brown PetscFunctionReturn(0); 118*c4762a1bSJed Brown } 119*c4762a1bSJed Brown 120*c4762a1bSJed Brown /* 121*c4762a1bSJed Brown Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian. 122*c4762a1bSJed Brown */ 123*c4762a1bSJed Brown PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx) 124*c4762a1bSJed Brown { 125*c4762a1bSJed Brown PetscErrorCode ierr; 126*c4762a1bSJed Brown PetscInt rowcol[] = {0,1,2}; 127*c4762a1bSJed Brown PetscScalar J[3][3]; 128*c4762a1bSJed Brown const PetscScalar *u,*udot; 129*c4762a1bSJed Brown 130*c4762a1bSJed Brown PetscFunctionBegin; 131*c4762a1bSJed Brown ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 132*c4762a1bSJed Brown ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 133*c4762a1bSJed Brown J[0][0] = a + ctx->k*u[1]; J[0][1] = ctx->k*u[0]; J[0][2] = 0.0; 134*c4762a1bSJed Brown J[1][0] = ctx->k*u[1]; J[1][1] = a + ctx->k*u[0]; J[1][2] = 0.0; 135*c4762a1bSJed Brown J[2][0] = -ctx->k*u[1]; J[2][1] = -ctx->k*u[0]; J[2][2] = a; 136*c4762a1bSJed Brown ierr = MatSetValues(B,3,rowcol,3,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 137*c4762a1bSJed Brown ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 138*c4762a1bSJed Brown ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 139*c4762a1bSJed Brown 140*c4762a1bSJed Brown ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 141*c4762a1bSJed Brown ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 142*c4762a1bSJed Brown if (A != B) { 143*c4762a1bSJed Brown ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 144*c4762a1bSJed Brown ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 145*c4762a1bSJed Brown } 146*c4762a1bSJed Brown PetscFunctionReturn(0); 147*c4762a1bSJed Brown } 148*c4762a1bSJed Brown 149*c4762a1bSJed Brown /* 150*c4762a1bSJed Brown Defines the exact (analytic) solution to the ODE 151*c4762a1bSJed Brown */ 152*c4762a1bSJed Brown static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx) 153*c4762a1bSJed Brown { 154*c4762a1bSJed Brown PetscErrorCode ierr; 155*c4762a1bSJed Brown const PetscScalar *uinit; 156*c4762a1bSJed Brown PetscScalar *u,d0,q; 157*c4762a1bSJed Brown 158*c4762a1bSJed Brown PetscFunctionBegin; 159*c4762a1bSJed Brown ierr = VecGetArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 160*c4762a1bSJed Brown ierr = VecGetArray(U,&u);CHKERRQ(ierr); 161*c4762a1bSJed Brown d0 = uinit[0] - uinit[1]; 162*c4762a1bSJed Brown if (d0 == 0.0) q = ctx->k*t; 163*c4762a1bSJed Brown else q = (1.0 - PetscExpScalar(-ctx->k*t*d0))/d0; 164*c4762a1bSJed Brown u[0] = uinit[0]/(1.0 + uinit[1]*q); 165*c4762a1bSJed Brown u[1] = u[0] - d0; 166*c4762a1bSJed Brown u[2] = uinit[1] + uinit[2] - u[1]; 167*c4762a1bSJed Brown ierr = VecRestoreArray(U,&u);CHKERRQ(ierr); 168*c4762a1bSJed Brown ierr = VecRestoreArrayRead(ctx->initialsolution,&uinit);CHKERRQ(ierr); 169*c4762a1bSJed Brown PetscFunctionReturn(0); 170*c4762a1bSJed Brown } 171*c4762a1bSJed Brown 172*c4762a1bSJed Brown int main(int argc,char **argv) 173*c4762a1bSJed Brown { 174*c4762a1bSJed Brown TS ts; /* ODE integrator */ 175*c4762a1bSJed Brown Vec U; /* solution will be stored here */ 176*c4762a1bSJed Brown Mat A; /* Jacobian matrix */ 177*c4762a1bSJed Brown PetscErrorCode ierr; 178*c4762a1bSJed Brown PetscMPIInt size; 179*c4762a1bSJed Brown PetscInt n = 3; 180*c4762a1bSJed Brown AppCtx ctx; 181*c4762a1bSJed Brown PetscScalar *u; 182*c4762a1bSJed Brown const char * const names[] = {"U1","U2","U3",NULL}; 183*c4762a1bSJed Brown 184*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 185*c4762a1bSJed Brown Initialize program 186*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 187*c4762a1bSJed Brown ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 188*c4762a1bSJed Brown ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRQ(ierr); 189*c4762a1bSJed Brown if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 190*c4762a1bSJed Brown 191*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 192*c4762a1bSJed Brown Create necessary matrix and vectors 193*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 194*c4762a1bSJed Brown ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 195*c4762a1bSJed Brown ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 196*c4762a1bSJed Brown ierr = MatSetFromOptions(A);CHKERRQ(ierr); 197*c4762a1bSJed Brown ierr = MatSetUp(A);CHKERRQ(ierr); 198*c4762a1bSJed Brown 199*c4762a1bSJed Brown ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 200*c4762a1bSJed Brown 201*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 202*c4762a1bSJed Brown Set runtime options 203*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 204*c4762a1bSJed Brown ctx.k = .9; 205*c4762a1bSJed Brown ierr = PetscOptionsGetScalar(NULL,NULL,"-k",&ctx.k,NULL);CHKERRQ(ierr); 206*c4762a1bSJed Brown ierr = VecDuplicate(U,&ctx.initialsolution);CHKERRQ(ierr); 207*c4762a1bSJed Brown ierr = VecGetArray(ctx.initialsolution,&u);CHKERRQ(ierr); 208*c4762a1bSJed Brown u[0] = 1; 209*c4762a1bSJed Brown u[1] = .7; 210*c4762a1bSJed Brown u[2] = 0; 211*c4762a1bSJed Brown ierr = VecRestoreArray(ctx.initialsolution,&u);CHKERRQ(ierr); 212*c4762a1bSJed Brown ierr = PetscOptionsGetVec(NULL,NULL,"-initial",ctx.initialsolution,NULL);CHKERRQ(ierr); 213*c4762a1bSJed Brown 214*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 215*c4762a1bSJed Brown Create timestepping solver context 216*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 217*c4762a1bSJed Brown ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 218*c4762a1bSJed Brown ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 219*c4762a1bSJed Brown ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); 220*c4762a1bSJed Brown ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx);CHKERRQ(ierr); 221*c4762a1bSJed Brown ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);CHKERRQ(ierr); 222*c4762a1bSJed Brown ierr = TSSetSolutionFunction(ts,(TSSolutionFunction)Solution,&ctx);CHKERRQ(ierr); 223*c4762a1bSJed Brown 224*c4762a1bSJed Brown { 225*c4762a1bSJed Brown DM dm; 226*c4762a1bSJed Brown void *ptr; 227*c4762a1bSJed Brown ierr = TSGetDM(ts,&dm);CHKERRQ(ierr); 228*c4762a1bSJed Brown ierr = PetscDLSym(NULL,"IFunctionView",&ptr);CHKERRQ(ierr); 229*c4762a1bSJed Brown ierr = PetscDLSym(NULL,"IFunctionLoad",&ptr);CHKERRQ(ierr); 230*c4762a1bSJed Brown ierr = DMTSSetIFunctionSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 231*c4762a1bSJed Brown ierr = DMTSSetIJacobianSerialize(dm,(PetscErrorCode (*)(void*,PetscViewer))IFunctionView,(PetscErrorCode (*)(void**,PetscViewer))IFunctionLoad);CHKERRQ(ierr); 232*c4762a1bSJed Brown } 233*c4762a1bSJed Brown 234*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 235*c4762a1bSJed Brown Set initial conditions 236*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 237*c4762a1bSJed Brown ierr = Solution(ts,0,U,&ctx);CHKERRQ(ierr); 238*c4762a1bSJed Brown ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 239*c4762a1bSJed Brown 240*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 241*c4762a1bSJed Brown Set solver options 242*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 243*c4762a1bSJed Brown ierr = TSSetTimeStep(ts,.001);CHKERRQ(ierr); 244*c4762a1bSJed Brown ierr = TSSetMaxSteps(ts,1000);CHKERRQ(ierr); 245*c4762a1bSJed Brown ierr = TSSetMaxTime(ts,20.0);CHKERRQ(ierr); 246*c4762a1bSJed Brown ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 247*c4762a1bSJed Brown ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 248*c4762a1bSJed Brown ierr = TSMonitorLGSetVariableNames(ts,names);CHKERRQ(ierr); 249*c4762a1bSJed Brown 250*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 251*c4762a1bSJed Brown Solve nonlinear system 252*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 253*c4762a1bSJed Brown ierr = TSSolve(ts,U);CHKERRQ(ierr); 254*c4762a1bSJed Brown 255*c4762a1bSJed Brown ierr = TSView(ts,PETSC_VIEWER_BINARY_WORLD);CHKERRQ(ierr); 256*c4762a1bSJed Brown 257*c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 258*c4762a1bSJed Brown Free work space. All PETSc objects should be destroyed when they are no longer needed. 259*c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 260*c4762a1bSJed Brown ierr = VecDestroy(&ctx.initialsolution);CHKERRQ(ierr); 261*c4762a1bSJed Brown ierr = MatDestroy(&A);CHKERRQ(ierr); 262*c4762a1bSJed Brown ierr = VecDestroy(&U);CHKERRQ(ierr); 263*c4762a1bSJed Brown ierr = TSDestroy(&ts);CHKERRQ(ierr); 264*c4762a1bSJed Brown 265*c4762a1bSJed Brown ierr = PetscFinalize(); 266*c4762a1bSJed Brown return ierr; 267*c4762a1bSJed Brown } 268*c4762a1bSJed Brown 269*c4762a1bSJed Brown 270*c4762a1bSJed Brown /*TEST 271*c4762a1bSJed Brown 272*c4762a1bSJed Brown test: 273*c4762a1bSJed Brown args: -ts_view 274*c4762a1bSJed Brown requires: dlsym define(PETSC_HAVE_DYNAMIC_LIBRARIES) 275*c4762a1bSJed Brown 276*c4762a1bSJed Brown test: 277*c4762a1bSJed Brown suffix: 2 278*c4762a1bSJed Brown args: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view 279*c4762a1bSJed Brown requires: x 280*c4762a1bSJed Brown output_file: output/ex1_1.out 281*c4762a1bSJed Brown requires: dlsym define(PETSC_HAVE_DYNAMIC_LIBRARIES) 282*c4762a1bSJed Brown 283*c4762a1bSJed Brown TEST*/ 284