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*/