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