Actual source code: ex2.c
2: static char help[] = "Reaction Equation from Chemistry\n";
4: /*
6: Page 6, An example from Atomospheric Chemistry
8: u_1_t =
9: u_2_t =
10: u_3_t =
11: u_4_t =
13: -ts_monitor_lg_error -ts_monitor_lg_solution -ts_view -ts_max_time 2.e4
15: */
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: */
27: #include <petscts.h>
29: typedef struct {
30: PetscScalar k1, k2, k3;
31: PetscScalar sigma2;
32: Vec initialsolution;
33: } AppCtx;
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: }
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;
48: VecGetArrayRead(U, &u);
49: VecGetArrayRead(Udot, &udot);
50: 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: VecRestoreArrayRead(U, &u);
56: VecRestoreArrayRead(Udot, &udot);
57: VecRestoreArrayWrite(F, &f);
58: return 0;
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: VecGetArrayRead(U, &u);
68: VecGetArrayRead(Udot, &udot);
69: J[0][0] = a + ctx->k2;
70: J[0][1] = 0.0;
71: J[0][2] = -k1(ctx, t);
72: J[0][3] = 0.0;
73: J[1][0] = 0.0;
74: J[1][1] = a + ctx->k3 * u[3];
75: J[1][2] = -k1(ctx, t);
76: J[1][3] = ctx->k3 * u[1];
77: J[2][0] = 0.0;
78: J[2][1] = -ctx->k3 * u[3];
79: J[2][2] = a + k1(ctx, t);
80: J[2][3] = -ctx->k3 * u[1];
81: J[3][0] = -ctx->k2;
82: J[3][1] = ctx->k3 * u[3];
83: J[3][2] = 0.0;
84: J[3][3] = a + ctx->k3 * u[1];
85: MatSetValues(B, 4, rowcol, 4, rowcol, &J[0][0], INSERT_VALUES);
86: VecRestoreArrayRead(U, &u);
87: VecRestoreArrayRead(Udot, &udot);
89: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
90: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
91: if (A != B) {
92: MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY);
93: MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY);
94: }
95: return 0;
96: }
98: static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx)
99: {
100: VecCopy(ctx->initialsolution, U);
102: return 0;
103: }
105: int main(int argc, char **argv)
106: {
107: TS ts; /* ODE integrator */
108: Vec U; /* solution */
109: Mat A; /* Jacobian matrix */
110: PetscMPIInt size;
111: PetscInt n = 4;
112: AppCtx ctx;
113: PetscScalar *u;
115: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
116: Initialize program
117: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
119: PetscInitialize(&argc, &argv, (char *)0, help);
120: MPI_Comm_size(PETSC_COMM_WORLD, &size);
123: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
124: Create necessary matrix and vectors
125: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
126: MatCreate(PETSC_COMM_WORLD, &A);
127: MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE);
128: MatSetFromOptions(A);
129: MatSetUp(A);
131: MatCreateVecs(A, &U, NULL);
133: ctx.k1 = 1.0e-5;
134: ctx.k2 = 1.0e5;
135: ctx.k3 = 1.0e-16;
136: ctx.sigma2 = 1.0e6;
138: VecDuplicate(U, &ctx.initialsolution);
139: VecGetArrayWrite(ctx.initialsolution, &u);
140: u[0] = 0.0;
141: u[1] = 1.3e8;
142: u[2] = 5.0e11;
143: u[3] = 8.0e11;
144: VecRestoreArrayWrite(ctx.initialsolution, &u);
146: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
147: Create timestepping solver context
148: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
149: TSCreate(PETSC_COMM_WORLD, &ts);
150: TSSetProblemType(ts, TS_NONLINEAR);
151: TSSetType(ts, TSROSW);
152: TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &ctx);
153: TSSetIJacobian(ts, A, A, (TSIJacobian)IJacobian, &ctx);
155: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
156: Set initial conditions
157: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
158: Solution(ts, 0, U, &ctx);
159: TSSetTime(ts, 4.0 * 3600);
160: TSSetTimeStep(ts, 1.0);
161: TSSetSolution(ts, U);
163: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
164: Set solver options
165: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
166: TSSetMaxTime(ts, 518400.0);
167: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
168: TSSetMaxStepRejections(ts, 100);
169: TSSetMaxSNESFailures(ts, -1); /* unlimited */
170: TSSetFromOptions(ts);
172: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
173: Solve nonlinear system
174: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
175: TSSolve(ts, U);
177: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
178: Free work space. All PETSc objects should be destroyed when they
179: are no longer needed.
180: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
181: VecDestroy(&ctx.initialsolution);
182: MatDestroy(&A);
183: VecDestroy(&U);
184: TSDestroy(&ts);
186: PetscFinalize();
187: return 0;
188: }
190: /*TEST
192: test:
193: args: -ts_view -ts_max_time 2.e4
194: timeoutfactor: 15
195: requires: !single
197: TEST*/