1 2 static char help[] = "Reaction Equation from Chemistry\n"; 3 4 /* 5 6 Page 6, An example from Atomospheric Chemistry 7 8 u_1_t = 9 u_2_t = 10 u_3_t = 11 u_4_t = 12 13 -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view -ts_max_time 2.e4 14 15 */ 16 17 /* 18 Include "petscts.h" so that we can use TS solvers. Note that this 19 file automatically includes: 20 petscsys.h - base PETSc routines petscvec.h - vectors 21 petscmat.h - matrices 22 petscis.h - index sets petscksp.h - Krylov subspace methods 23 petscviewer.h - viewers petscpc.h - preconditioners 24 petscksp.h - linear solvers 25 */ 26 27 #include <petscts.h> 28 29 typedef struct { 30 PetscScalar k1,k2,k3; 31 PetscScalar sigma2; 32 Vec initialsolution; 33 } AppCtx; 34 35 PetscScalar k1(AppCtx *ctx,PetscReal t) 36 { 37 PetscReal th = t/3600.0; 38 PetscReal barth = th - 24.0*PetscFloorReal(th/24.0); 39 if (((((PetscInt)th) % 24) < 4) || ((((PetscInt)th) % 24) >= 20)) return(1.0e-40); 40 else return(ctx->k1*PetscExpReal(7.0*PetscPowReal(PetscSinReal(.0625*PETSC_PI*(barth - 4.0)),.2))); 41 } 42 43 static PetscErrorCode IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx *ctx) 44 { 45 PetscErrorCode ierr; 46 PetscScalar *f; 47 const PetscScalar *u,*udot; 48 49 PetscFunctionBegin; 50 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 51 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 52 ierr = VecGetArrayWrite(F,&f);CHKERRQ(ierr); 53 f[0] = udot[0] - k1(ctx,t)*u[2] + ctx->k2*u[0]; 54 f[1] = udot[1] - k1(ctx,t)*u[2] + ctx->k3*u[1]*u[3] - ctx->sigma2; 55 f[2] = udot[2] - ctx->k3*u[1]*u[3] + k1(ctx,t)*u[2]; 56 f[3] = udot[3] - ctx->k2*u[0] + ctx->k3*u[1]*u[3]; 57 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 58 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 59 ierr = VecRestoreArrayWrite(F,&f);CHKERRQ(ierr); 60 PetscFunctionReturn(0); 61 } 62 63 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx) 64 { 65 PetscErrorCode ierr; 66 PetscInt rowcol[] = {0,1,2,3}; 67 PetscScalar J[4][4]; 68 const PetscScalar *u,*udot; 69 70 PetscFunctionBegin; 71 ierr = VecGetArrayRead(U,&u);CHKERRQ(ierr); 72 ierr = VecGetArrayRead(Udot,&udot);CHKERRQ(ierr); 73 J[0][0] = a + ctx->k2; J[0][1] = 0.0; J[0][2] = -k1(ctx,t); J[0][3] = 0.0; 74 J[1][0] = 0.0; J[1][1] = a + ctx->k3*u[3]; J[1][2] = -k1(ctx,t); J[1][3] = ctx->k3*u[1]; 75 J[2][0] = 0.0; J[2][1] = -ctx->k3*u[3]; J[2][2] = a + k1(ctx,t); J[2][3] = -ctx->k3*u[1]; 76 J[3][0] = -ctx->k2; J[3][1] = ctx->k3*u[3]; J[3][2] = 0.0; J[3][3] = a + ctx->k3*u[1]; 77 ierr = MatSetValues(B,4,rowcol,4,rowcol,&J[0][0],INSERT_VALUES);CHKERRQ(ierr); 78 ierr = VecRestoreArrayRead(U,&u);CHKERRQ(ierr); 79 ierr = VecRestoreArrayRead(Udot,&udot);CHKERRQ(ierr); 80 81 ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 82 ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 83 if (A != B) { 84 ierr = MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 85 ierr = MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY);CHKERRQ(ierr); 86 } 87 PetscFunctionReturn(0); 88 } 89 90 static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx) 91 { 92 PetscErrorCode ierr; 93 94 PetscFunctionBegin; 95 ierr = VecCopy(ctx->initialsolution,U);CHKERRQ(ierr); 96 if (t > 0) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Solution not given"); 97 PetscFunctionReturn(0); 98 } 99 100 int main(int argc,char **argv) 101 { 102 TS ts; /* ODE integrator */ 103 Vec U; /* solution */ 104 Mat A; /* Jacobian matrix */ 105 PetscErrorCode ierr; 106 PetscMPIInt size; 107 PetscInt n = 4; 108 AppCtx ctx; 109 PetscScalar *u; 110 111 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 112 Initialize program 113 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 114 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr; 115 ierr = MPI_Comm_size(PETSC_COMM_WORLD,&size);CHKERRMPI(ierr); 116 if (size > 1) SETERRQ(PETSC_COMM_WORLD,PETSC_ERR_SUP,"Only for sequential runs"); 117 118 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 119 Create necessary matrix and vectors 120 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 121 ierr = MatCreate(PETSC_COMM_WORLD,&A);CHKERRQ(ierr); 122 ierr = MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE);CHKERRQ(ierr); 123 ierr = MatSetFromOptions(A);CHKERRQ(ierr); 124 ierr = MatSetUp(A);CHKERRQ(ierr); 125 126 ierr = MatCreateVecs(A,&U,NULL);CHKERRQ(ierr); 127 128 ctx.k1 = 1.0e-5; 129 ctx.k2 = 1.0e5; 130 ctx.k3 = 1.0e-16; 131 ctx.sigma2 = 1.0e6; 132 133 ierr = VecDuplicate(U,&ctx.initialsolution);CHKERRQ(ierr); 134 ierr = VecGetArrayWrite(ctx.initialsolution,&u);CHKERRQ(ierr); 135 u[0] = 0.0; 136 u[1] = 1.3e8; 137 u[2] = 5.0e11; 138 u[3] = 8.0e11; 139 ierr = VecRestoreArrayWrite(ctx.initialsolution,&u);CHKERRQ(ierr); 140 141 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 142 Create timestepping solver context 143 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 144 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr); 145 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr); 146 ierr = TSSetType(ts,TSROSW);CHKERRQ(ierr); 147 ierr = TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx);CHKERRQ(ierr); 148 ierr = TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx);CHKERRQ(ierr); 149 150 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 151 Set initial conditions 152 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 153 ierr = Solution(ts,0,U,&ctx);CHKERRQ(ierr); 154 ierr = TSSetTime(ts,4.0*3600);CHKERRQ(ierr); 155 ierr = TSSetTimeStep(ts,1.0);CHKERRQ(ierr); 156 ierr = TSSetSolution(ts,U);CHKERRQ(ierr); 157 158 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 159 Set solver options 160 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 161 ierr = TSSetMaxTime(ts,518400.0);CHKERRQ(ierr); 162 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr); 163 ierr = TSSetMaxStepRejections(ts,100);CHKERRQ(ierr); 164 ierr = TSSetMaxSNESFailures(ts,-1);CHKERRQ(ierr); /* unlimited */ 165 ierr = TSSetFromOptions(ts);CHKERRQ(ierr); 166 167 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 168 Solve nonlinear system 169 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 170 ierr = TSSolve(ts,U);CHKERRQ(ierr); 171 172 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 173 Free work space. All PETSc objects should be destroyed when they 174 are no longer needed. 175 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 176 ierr = VecDestroy(&ctx.initialsolution);CHKERRQ(ierr); 177 ierr = MatDestroy(&A);CHKERRQ(ierr); 178 ierr = VecDestroy(&U);CHKERRQ(ierr); 179 ierr = TSDestroy(&ts);CHKERRQ(ierr); 180 181 ierr = PetscFinalize(); 182 return ierr; 183 } 184 185 /*TEST 186 187 test: 188 args: -ts_view -ts_max_time 2.e4 189 timeoutfactor: 15 190 requires: !single 191 192 TEST*/ 193