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 PetscScalar *f; 46 const PetscScalar *u,*udot; 47 48 PetscFunctionBegin; 49 PetscCall(VecGetArrayRead(U,&u)); 50 PetscCall(VecGetArrayRead(Udot,&udot)); 51 PetscCall(VecGetArrayWrite(F,&f)); 52 f[0] = udot[0] - k1(ctx,t)*u[2] + ctx->k2*u[0]; 53 f[1] = udot[1] - k1(ctx,t)*u[2] + ctx->k3*u[1]*u[3] - ctx->sigma2; 54 f[2] = udot[2] - ctx->k3*u[1]*u[3] + k1(ctx,t)*u[2]; 55 f[3] = udot[3] - ctx->k2*u[0] + ctx->k3*u[1]*u[3]; 56 PetscCall(VecRestoreArrayRead(U,&u)); 57 PetscCall(VecRestoreArrayRead(Udot,&udot)); 58 PetscCall(VecRestoreArrayWrite(F,&f)); 59 PetscFunctionReturn(0); 60 } 61 62 static PetscErrorCode IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx *ctx) 63 { 64 PetscInt rowcol[] = {0,1,2,3}; 65 PetscScalar J[4][4]; 66 const PetscScalar *u,*udot; 67 68 PetscFunctionBegin; 69 PetscCall(VecGetArrayRead(U,&u)); 70 PetscCall(VecGetArrayRead(Udot,&udot)); 71 J[0][0] = a + ctx->k2; J[0][1] = 0.0; J[0][2] = -k1(ctx,t); J[0][3] = 0.0; 72 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]; 73 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]; 74 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]; 75 PetscCall(MatSetValues(B,4,rowcol,4,rowcol,&J[0][0],INSERT_VALUES)); 76 PetscCall(VecRestoreArrayRead(U,&u)); 77 PetscCall(VecRestoreArrayRead(Udot,&udot)); 78 79 PetscCall(MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY)); 80 PetscCall(MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY)); 81 if (A != B) { 82 PetscCall(MatAssemblyBegin(B,MAT_FINAL_ASSEMBLY)); 83 PetscCall(MatAssemblyEnd(B,MAT_FINAL_ASSEMBLY)); 84 } 85 PetscFunctionReturn(0); 86 } 87 88 static PetscErrorCode Solution(TS ts,PetscReal t,Vec U,AppCtx *ctx) 89 { 90 PetscFunctionBegin; 91 PetscCall(VecCopy(ctx->initialsolution,U)); 92 PetscCheck(t <= 0,PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Solution not given"); 93 PetscFunctionReturn(0); 94 } 95 96 int main(int argc,char **argv) 97 { 98 TS ts; /* ODE integrator */ 99 Vec U; /* solution */ 100 Mat A; /* Jacobian matrix */ 101 PetscMPIInt size; 102 PetscInt n = 4; 103 AppCtx ctx; 104 PetscScalar *u; 105 106 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 107 Initialize program 108 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 109 PetscFunctionBeginUser; 110 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 111 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 112 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 113 114 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 115 Create necessary matrix and vectors 116 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 117 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 118 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 119 PetscCall(MatSetFromOptions(A)); 120 PetscCall(MatSetUp(A)); 121 122 PetscCall(MatCreateVecs(A,&U,NULL)); 123 124 ctx.k1 = 1.0e-5; 125 ctx.k2 = 1.0e5; 126 ctx.k3 = 1.0e-16; 127 ctx.sigma2 = 1.0e6; 128 129 PetscCall(VecDuplicate(U,&ctx.initialsolution)); 130 PetscCall(VecGetArrayWrite(ctx.initialsolution,&u)); 131 u[0] = 0.0; 132 u[1] = 1.3e8; 133 u[2] = 5.0e11; 134 u[3] = 8.0e11; 135 PetscCall(VecRestoreArrayWrite(ctx.initialsolution,&u)); 136 137 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 138 Create timestepping solver context 139 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 140 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 141 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 142 PetscCall(TSSetType(ts,TSROSW)); 143 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 144 PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 145 146 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 147 Set initial conditions 148 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 149 PetscCall(Solution(ts,0,U,&ctx)); 150 PetscCall(TSSetTime(ts,4.0*3600)); 151 PetscCall(TSSetTimeStep(ts,1.0)); 152 PetscCall(TSSetSolution(ts,U)); 153 154 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 155 Set solver options 156 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 157 PetscCall(TSSetMaxTime(ts,518400.0)); 158 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 159 PetscCall(TSSetMaxStepRejections(ts,100)); 160 PetscCall(TSSetMaxSNESFailures(ts,-1)); /* unlimited */ 161 PetscCall(TSSetFromOptions(ts)); 162 163 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 164 Solve nonlinear system 165 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 166 PetscCall(TSSolve(ts,U)); 167 168 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 169 Free work space. All PETSc objects should be destroyed when they 170 are no longer needed. 171 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 172 PetscCall(VecDestroy(&ctx.initialsolution)); 173 PetscCall(MatDestroy(&A)); 174 PetscCall(VecDestroy(&U)); 175 PetscCall(TSDestroy(&ts)); 176 177 PetscCall(PetscFinalize()); 178 return 0; 179 } 180 181 /*TEST 182 183 test: 184 args: -ts_view -ts_max_time 2.e4 185 timeoutfactor: 15 186 requires: !single 187 188 TEST*/ 189