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:   PetscFunctionBegin;
 49:   PetscCall(VecGetArrayRead(U, &u));
 50:   PetscCall(VecGetArrayRead(Udot, &udot));
 51:   PetscCall(VecGetArrayWrite(F, &f));
 52:   f[0] = udot[0] - k1(ctx, t) * u[2] + ctx->k2 * u[0];
 53:   f[1] = udot[1] - k1(ctx, t) * u[2] + ctx->k3 * u[1] * u[3] - ctx->sigma2;
 54:   f[2] = udot[2] - ctx->k3 * u[1] * u[3] + k1(ctx, t) * u[2];
 55:   f[3] = udot[3] - ctx->k2 * u[0] + ctx->k3 * u[1] * u[3];
 56:   PetscCall(VecRestoreArrayRead(U, &u));
 57:   PetscCall(VecRestoreArrayRead(Udot, &udot));
 58:   PetscCall(VecRestoreArrayWrite(F, &f));
 59:   PetscFunctionReturn(PETSC_SUCCESS);
 60: }

 62: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, AppCtx *ctx)
 63: {
 64:   PetscInt           rowcol[] = {0, 1, 2, 3};
 65:   PetscScalar        J[4][4];
 66:   const PetscScalar *u, *udot;

 68:   PetscFunctionBegin;
 69:   PetscCall(VecGetArrayRead(U, &u));
 70:   PetscCall(VecGetArrayRead(Udot, &udot));
 71:   J[0][0] = a + ctx->k2;
 72:   J[0][1] = 0.0;
 73:   J[0][2] = -k1(ctx, t);
 74:   J[0][3] = 0.0;
 75:   J[1][0] = 0.0;
 76:   J[1][1] = a + ctx->k3 * u[3];
 77:   J[1][2] = -k1(ctx, t);
 78:   J[1][3] = ctx->k3 * u[1];
 79:   J[2][0] = 0.0;
 80:   J[2][1] = -ctx->k3 * u[3];
 81:   J[2][2] = a + k1(ctx, t);
 82:   J[2][3] = -ctx->k3 * u[1];
 83:   J[3][0] = -ctx->k2;
 84:   J[3][1] = ctx->k3 * u[3];
 85:   J[3][2] = 0.0;
 86:   J[3][3] = a + ctx->k3 * u[1];
 87:   PetscCall(MatSetValues(B, 4, rowcol, 4, rowcol, &J[0][0], INSERT_VALUES));
 88:   PetscCall(VecRestoreArrayRead(U, &u));
 89:   PetscCall(VecRestoreArrayRead(Udot, &udot));

 91:   PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
 92:   PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
 93:   if (A != B) {
 94:     PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
 95:     PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
 96:   }
 97:   PetscFunctionReturn(PETSC_SUCCESS);
 98: }

100: static PetscErrorCode Solution(TS ts, PetscReal t, Vec U, AppCtx *ctx)
101: {
102:   PetscFunctionBegin;
103:   PetscCall(VecCopy(ctx->initialsolution, U));
104:   PetscCheck(t <= 0, PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "Solution not given");
105:   PetscFunctionReturn(PETSC_SUCCESS);
106: }

108: int main(int argc, char **argv)
109: {
110:   TS           ts; /* ODE integrator */
111:   Vec          U;  /* solution */
112:   Mat          A;  /* Jacobian matrix */
113:   PetscMPIInt  size;
114:   PetscInt     n = 4;
115:   AppCtx       ctx;
116:   PetscScalar *u;

118:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
119:      Initialize program
120:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
121:   PetscFunctionBeginUser;
122:   PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
123:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
124:   PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");

126:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
127:     Create necessary matrix and vectors
128:     - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
129:   PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
130:   PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
131:   PetscCall(MatSetFromOptions(A));
132:   PetscCall(MatSetUp(A));

134:   PetscCall(MatCreateVecs(A, &U, NULL));

136:   ctx.k1     = 1.0e-5;
137:   ctx.k2     = 1.0e5;
138:   ctx.k3     = 1.0e-16;
139:   ctx.sigma2 = 1.0e6;

141:   PetscCall(VecDuplicate(U, &ctx.initialsolution));
142:   PetscCall(VecGetArrayWrite(ctx.initialsolution, &u));
143:   u[0] = 0.0;
144:   u[1] = 1.3e8;
145:   u[2] = 5.0e11;
146:   u[3] = 8.0e11;
147:   PetscCall(VecRestoreArrayWrite(ctx.initialsolution, &u));

149:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
150:      Create timestepping solver context
151:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
152:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
153:   PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
154:   PetscCall(TSSetType(ts, TSROSW));
155:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &ctx));
156:   PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobian)IJacobian, &ctx));

158:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
159:      Set initial conditions
160:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
161:   PetscCall(Solution(ts, 0, U, &ctx));
162:   PetscCall(TSSetTime(ts, 4.0 * 3600));
163:   PetscCall(TSSetTimeStep(ts, 1.0));
164:   PetscCall(TSSetSolution(ts, U));

166:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167:      Set solver options
168:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169:   PetscCall(TSSetMaxTime(ts, 518400.0));
170:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
171:   PetscCall(TSSetMaxStepRejections(ts, 100));
172:   PetscCall(TSSetMaxSNESFailures(ts, -1)); /* unlimited */
173:   PetscCall(TSSetFromOptions(ts));

175:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
176:      Solve nonlinear system
177:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
178:   PetscCall(TSSolve(ts, U));

180:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
181:      Free work space.  All PETSc objects should be destroyed when they
182:      are no longer needed.
183:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
184:   PetscCall(VecDestroy(&ctx.initialsolution));
185:   PetscCall(MatDestroy(&A));
186:   PetscCall(VecDestroy(&U));
187:   PetscCall(TSDestroy(&ts));

189:   PetscCall(PetscFinalize());
190:   return 0;
191: }

193: /*TEST

195:    test:
196:      args: -ts_view -ts_max_time 2.e4
197:      timeoutfactor: 15
198:      requires: !single

200: TEST*/