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