Actual source code: ex51.c

  1: static char help[] = "Small ODE to test TS accuracy.\n";

  3: /*
  4:   The ODE
  5:                   u1_t = cos(t),
  6:                   u2_t = sin(u2)
  7:   with analytical solution
  8:                   u1(t) = sin(t),
  9:                   u2(t) = 2 * atan(exp(t) * tan(0.5))
 10:   is used to test the accuracy of TS schemes.
 11: */

 13: #include <petscts.h>

 15: /*
 16:      Defines the ODE passed to the ODE solver in explicit form: U_t = F(U)
 17: */
 18: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, void *s)
 19: {
 20:   PetscScalar       *f;
 21:   const PetscScalar *u;

 24:   VecGetArrayRead(U, &u);
 25:   VecGetArray(F, &f);

 27:   f[0] = PetscCosReal(t);
 28:   f[1] = PetscSinReal(u[1]);

 30:   VecRestoreArrayRead(U, &u);
 31:   VecRestoreArray(F, &f);
 32:   return 0;
 33: }

 35: /*
 36:      Defines the exact solution.
 37: */
 38: static PetscErrorCode ExactSolution(PetscReal t, Vec U)
 39: {
 40:   PetscScalar *u;

 43:   VecGetArray(U, &u);

 45:   u[0] = PetscSinReal(t);
 46:   u[1] = 2 * PetscAtanReal(PetscExpReal(t) * PetscTanReal(0.5));

 48:   VecRestoreArray(U, &u);
 49:   return 0;
 50: }

 52: int main(int argc, char **argv)
 53: {
 54:   TS           ts;  /* ODE integrator */
 55:   Vec          U;   /* numerical solution will be stored here */
 56:   Vec          Uex; /* analytical (exact) solution will be stored here */
 57:   PetscMPIInt  size;
 58:   PetscInt     n = 2;
 59:   PetscScalar *u;
 60:   PetscReal    t, final_time = 1.0, dt = 0.25;
 61:   PetscReal    error;
 62:   TSAdapt      adapt;

 64:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 65:      Initialize program
 66:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 68:   PetscInitialize(&argc, &argv, (char *)0, help);
 69:   MPI_Comm_size(PETSC_COMM_WORLD, &size);

 72:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 73:      Create timestepping solver context
 74:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 75:   TSCreate(PETSC_COMM_WORLD, &ts);
 76:   TSSetType(ts, TSROSW);

 78:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 79:      Set ODE routines
 80:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 81:   TSSetProblemType(ts, TS_NONLINEAR);
 82:   TSSetRHSFunction(ts, NULL, RHSFunction, NULL);

 84:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 85:      Set initial conditions
 86:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 87:   VecCreate(PETSC_COMM_WORLD, &U);
 88:   VecSetSizes(U, n, PETSC_DETERMINE);
 89:   VecSetUp(U);
 90:   VecGetArray(U, &u);
 91:   u[0] = 0.0;
 92:   u[1] = 1.0;
 93:   VecRestoreArray(U, &u);
 94:   TSSetSolution(ts, U);

 96:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 97:      Set solver options
 98:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
 99:   TSSetSaveTrajectory(ts);
100:   TSSetMaxTime(ts, final_time);
101:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
102:   TSSetTimeStep(ts, dt);
103:   /* The adaptive time step controller is forced to take constant time steps. */
104:   TSGetAdapt(ts, &adapt);
105:   TSAdaptSetType(adapt, TSADAPTNONE);

107:   TSSetFromOptions(ts);

109:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
110:      Run timestepping solver and compute error
111:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
112:   TSSolve(ts, U);
113:   TSGetTime(ts, &t);

115:   if (PetscAbsReal(t - final_time) > 100 * PETSC_MACHINE_EPSILON) {
116:     PetscPrintf(PETSC_COMM_WORLD, "Note: There is a difference of %g between the prescribed final time %g and the actual final time.\n", (double)(final_time - t), (double)final_time);
117:   }
118:   VecDuplicate(U, &Uex);
119:   ExactSolution(t, Uex);

121:   VecAYPX(Uex, -1.0, U);
122:   VecNorm(Uex, NORM_2, &error);
123:   PetscPrintf(PETSC_COMM_WORLD, "Error at final time: %.2E\n", (double)error);

125:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
126:      Free work space.  All PETSc objects should be destroyed when they are no longer needed.
127:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
128:   VecDestroy(&U);
129:   VecDestroy(&Uex);
130:   TSDestroy(&ts);

132:   PetscFinalize();
133:   return 0;
134: }

136: /*TEST

138:     test:
139:       suffix: 3bs
140:       args: -ts_type rk -ts_rk_type 3bs
141:       requires: !single

143:     test:
144:       suffix: 5bs
145:       args: -ts_type rk -ts_rk_type 5bs
146:       requires: !single

148:     test:
149:       suffix: 5dp
150:       args: -ts_type rk -ts_rk_type 5dp
151:       requires: !single

153:     test:
154:       suffix: 6vr
155:       args: -ts_type rk -ts_rk_type 6vr
156:       requires: !single

158:     test:
159:       suffix: 7vr
160:       args: -ts_type rk -ts_rk_type 7vr
161:       requires: !single

163:     test:
164:       suffix: 8vr
165:       args: -ts_type rk -ts_rk_type 8vr
166:       requires: !single

168: TEST*/