Actual source code: ex14.c

  1: static char help[] = "Tests all TSRK types \n\n";

  3: #include <petscts.h>

  5: static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ctx)
  6: {
  7:   PetscInt           i, n;
  8:   const PetscScalar *xx;
  9:   /* */ PetscScalar *ff;

 12:   VecGetLocalSize(X, &n);
 13:   VecGetArrayRead(X, &xx);
 14:   VecGetArray(F, &ff);

 16:   if (n >= 1) ff[0] = 1;
 17:   for (i = 1; i < n; i++) ff[i] = (i + 1) * (xx[i - 1] + PetscPowReal(t, i)) / 2;

 19:   VecRestoreArrayRead(X, &xx);
 20:   VecRestoreArray(F, &ff);
 21:   return 0;
 22: }

 24: PetscErrorCode TestCheckStage(TSAdapt adapt, TS ts, PetscReal t, Vec X, PetscBool *accept)
 25: {
 26:   PetscInt step;

 29:   TSGetStepNumber(ts, &step);
 30:   *accept = (step >= 2) ? PETSC_FALSE : PETSC_TRUE;
 31:   return 0;
 32: }

 34: static PetscErrorCode TestExplicitTS(TS ts, PetscInt order, const char subtype[])
 35: {
 36:   PetscInt           i;
 37:   PetscReal          t;
 38:   Vec                U, X, Y;
 39:   TSType             type;
 40:   PetscBool          done;
 41:   const PetscScalar *xx;
 42:   const PetscScalar *yy;
 43:   const PetscReal    Tf  = 1;
 44:   const PetscReal    dt  = Tf / 8;
 45:   const PetscReal    eps = 100 * PETSC_MACHINE_EPSILON;
 46:   TSAdapt            adapt;
 47:   PetscInt           step;
 48:   TSConvergedReason  reason;

 51:   TSGetType(ts, &type);
 52:   TSGetSolution(ts, &U);
 53:   VecZeroEntries(U);
 54:   TSSetStepNumber(ts, 0);
 55:   TSSetTime(ts, 0);
 56:   TSSetTimeStep(ts, dt);
 57:   TSSetMaxTime(ts, Tf);
 58:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
 59:   TSSolve(ts, NULL);
 60:   TSRollBack(ts);
 61:   TSSolve(ts, NULL);
 62:   TSGetTime(ts, &t);

 64:   TSGetSolution(ts, &U);
 65:   VecDuplicate(U, &X);
 66:   TSEvaluateStep(ts, order, X, NULL);
 67:   VecGetArrayRead(X, &xx);
 68:   for (i = 0; i < order; i++) {
 69:     PetscReal error = PetscAbsReal(PetscRealPart(xx[i]) - PetscPowReal(t, i + 1));
 71:   }
 72:   VecRestoreArrayRead(X, &xx);
 73:   VecDestroy(&X);

 75:   TSGetSolution(ts, &U);
 76:   VecDuplicate(U, &Y);
 77:   TSEvaluateStep(ts, order - 1, Y, &done);
 78:   VecGetArrayRead(Y, &yy);
 79:   for (i = 0; done && i < order - 1; i++) {
 80:     PetscReal error = PetscAbsReal(PetscRealPart(yy[i]) - PetscPowReal(t, i + 1));
 82:   }
 83:   VecRestoreArrayRead(Y, &yy);
 84:   VecDestroy(&Y);

 86:   TSGetAdapt(ts, &adapt);
 87:   TSAdaptSetCheckStage(adapt, TestCheckStage);
 88:   TSSetErrorIfStepFails(ts, PETSC_FALSE);
 89:   TSSetStepNumber(ts, 0);
 90:   TSSetTime(ts, 0);
 91:   TSSetTimeStep(ts, dt);
 92:   TSSolve(ts, NULL);
 93:   TSAdaptSetCheckStage(adapt, NULL);
 94:   TSSetErrorIfStepFails(ts, PETSC_TRUE);
 95:   TSGetStepNumber(ts, &step);
 96:   TSGetConvergedReason(ts, &reason);
 99:   return 0;
100: }

102: static PetscErrorCode TestTSRK(TS ts, TSRKType type)
103: {
104:   PetscInt    order;
105:   TSAdapt     adapt;
106:   PetscBool   rk1, rk3, rk4;
107:   TSAdaptType adapttype;
108:   char        savetype[32];

111:   TSRKSetType(ts, type);
112:   TSRKGetType(ts, &type);
113:   TSRKGetOrder(ts, &order);

115:   TSGetAdapt(ts, &adapt);
116:   TSAdaptGetType(adapt, &adapttype);
117:   PetscStrncpy(savetype, adapttype, sizeof(savetype));
118:   PetscStrcmp(type, TSRK1FE, &rk1);
119:   PetscStrcmp(type, TSRK3, &rk3);
120:   PetscStrcmp(type, TSRK4, &rk4);
121:   if (rk1 || rk3 || rk4) TSAdaptSetType(adapt, TSADAPTNONE);

123:   TestExplicitTS(ts, order, type);

125:   TSGetAdapt(ts, &adapt);
126:   TSAdaptSetType(adapt, savetype);
127:   return 0;
128: }

130: int main(int argc, char *argv[])
131: {
132:   TS       ts;
133:   Vec      X;
134:   PetscInt N = 9;

137:   PetscInitialize(&argc, &argv, NULL, help);

139:   TSCreate(PETSC_COMM_SELF, &ts);
140:   TSSetType(ts, TSRK);
141:   TSSetRHSFunction(ts, NULL, RHSFunction, NULL);
142:   VecCreateSeq(PETSC_COMM_SELF, N, &X);
143:   TSSetSolution(ts, X);
144:   VecDestroy(&X);
145:   TSSetFromOptions(ts);

147:   TestTSRK(ts, TSRK1FE);
148:   TestTSRK(ts, TSRK2A);
149:   TestTSRK(ts, TSRK3);
150:   TestTSRK(ts, TSRK3BS);
151:   TestTSRK(ts, TSRK4);
152:   TestTSRK(ts, TSRK5F);
153:   TestTSRK(ts, TSRK5DP);
154:   TestTSRK(ts, TSRK5BS);
155:   TestTSRK(ts, TSRK6VR);
156:   TestTSRK(ts, TSRK7VR);
157:   TestTSRK(ts, TSRK8VR);

159:   TSRollBack(ts);
160:   TSDestroy(&ts);

162:   PetscFinalize();
163:   return 0;
164: }

166: /*TEST

168: testset:
169:   output_file: output/ex14.out
170:   test:
171:     suffix: 0
172:   test:
173:     suffix: 1
174:     args: -ts_adapt_type none
175:   test:
176:     suffix: 2
177:     args: -ts_adapt_type basic
178:   test:
179:     suffix: 3
180:     args: -ts_adapt_type dsp

182: TEST*/