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