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
k1(AppCtx * ctx,PetscReal t)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
IFunction(TS ts,PetscReal t,Vec U,Vec Udot,Vec F,AppCtx * ctx)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
IJacobian(TS ts,PetscReal t,Vec U,Vec Udot,PetscReal a,Mat A,Mat B,AppCtx * ctx)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
Solution(TS ts,PetscReal t,Vec U,AppCtx * ctx)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
main(int argc,char ** argv)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, NULL, 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