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