Actual source code: ex2.c
1: static char help[] = "Reaction Equation from Chemistry\n";
3: /*
5: Page 6, An example from Atomospheric Chemistry
7: u_1_t =
8: u_2_t =
9: u_3_t =
10: u_4_t =
12: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view -ts_max_time 2.e4
14: */
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: */
26: #include <petscts.h>
28: typedef struct {
29: PetscScalar k1, k2, k3;
30: PetscScalar sigma2;
31: Vec initialsolution;
32: } AppCtx;
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: }
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;
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: }
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;
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));
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: }
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: }
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;
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");
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));
133: PetscCall(MatCreateVecs(A, &U, NULL));
135: ctx.k1 = 1.0e-5;
136: ctx.k2 = 1.0e5;
137: ctx.k3 = 1.0e-16;
138: ctx.sigma2 = 1.0e6;
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));
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));
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));
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));
174: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
175: Solve nonlinear system
176: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
177: PetscCall(TSSolve(ts, U));
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));
188: PetscCall(PetscFinalize());
189: return 0;
190: }
192: /*TEST
194: test:
195: args: -ts_view -ts_max_time 2.e4
196: timeoutfactor: 15
197: requires: !single
199: TEST*/