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 PetscCall(PetscInitialize(&argc,&argv,(char*)0,help)); 110 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 111 PetscCheck(size == 1,PETSC_COMM_WORLD,PETSC_ERR_WRONG_MPI_SIZE,"Only for sequential runs"); 112 113 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 114 Create necessary matrix and vectors 115 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 116 PetscCall(MatCreate(PETSC_COMM_WORLD,&A)); 117 PetscCall(MatSetSizes(A,n,n,PETSC_DETERMINE,PETSC_DETERMINE)); 118 PetscCall(MatSetFromOptions(A)); 119 PetscCall(MatSetUp(A)); 120 121 PetscCall(MatCreateVecs(A,&U,NULL)); 122 123 ctx.k1 = 1.0e-5; 124 ctx.k2 = 1.0e5; 125 ctx.k3 = 1.0e-16; 126 ctx.sigma2 = 1.0e6; 127 128 PetscCall(VecDuplicate(U,&ctx.initialsolution)); 129 PetscCall(VecGetArrayWrite(ctx.initialsolution,&u)); 130 u[0] = 0.0; 131 u[1] = 1.3e8; 132 u[2] = 5.0e11; 133 u[3] = 8.0e11; 134 PetscCall(VecRestoreArrayWrite(ctx.initialsolution,&u)); 135 136 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 137 Create timestepping solver context 138 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 139 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 140 PetscCall(TSSetProblemType(ts,TS_NONLINEAR)); 141 PetscCall(TSSetType(ts,TSROSW)); 142 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) IFunction,&ctx)); 143 PetscCall(TSSetIJacobian(ts,A,A,(TSIJacobian)IJacobian,&ctx)); 144 145 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 146 Set initial conditions 147 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 148 PetscCall(Solution(ts,0,U,&ctx)); 149 PetscCall(TSSetTime(ts,4.0*3600)); 150 PetscCall(TSSetTimeStep(ts,1.0)); 151 PetscCall(TSSetSolution(ts,U)); 152 153 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 154 Set solver options 155 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 156 PetscCall(TSSetMaxTime(ts,518400.0)); 157 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 158 PetscCall(TSSetMaxStepRejections(ts,100)); 159 PetscCall(TSSetMaxSNESFailures(ts,-1)); /* unlimited */ 160 PetscCall(TSSetFromOptions(ts)); 161 162 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 163 Solve nonlinear system 164 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 165 PetscCall(TSSolve(ts,U)); 166 167 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 168 Free work space. All PETSc objects should be destroyed when they 169 are no longer needed. 170 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 171 PetscCall(VecDestroy(&ctx.initialsolution)); 172 PetscCall(MatDestroy(&A)); 173 PetscCall(VecDestroy(&U)); 174 PetscCall(TSDestroy(&ts)); 175 176 PetscCall(PetscFinalize()); 177 return 0; 178 } 179 180 /*TEST 181 182 test: 183 args: -ts_view -ts_max_time 2.e4 184 timeoutfactor: 15 185 requires: !single 186 187 TEST*/ 188