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, PetscCtx ctx)
6: {
7: PetscInt n;
8: const PetscScalar *xx;
9: /* */ PetscScalar *ff;
11: PetscFunctionBeginUser;
12: PetscCall(VecGetLocalSize(X, &n));
13: PetscCall(VecGetArrayRead(X, &xx));
14: PetscCall(VecGetArray(F, &ff));
16: if (n >= 1) ff[0] = 1;
17: for (PetscInt i = 1; i < n; i++) ff[i] = (i + 1) * (xx[i - 1] + PetscPowReal(t, i)) / 2;
19: PetscCall(VecRestoreArrayRead(X, &xx));
20: PetscCall(VecRestoreArray(F, &ff));
21: PetscFunctionReturn(PETSC_SUCCESS);
22: }
24: PetscErrorCode TestCheckStage(TSAdapt adapt, TS ts, PetscReal t, Vec X, PetscBool *accept)
25: {
26: PetscInt step;
28: PetscFunctionBeginUser;
29: PetscCall(TSGetStepNumber(ts, &step));
30: *accept = (step >= 2) ? PETSC_FALSE : PETSC_TRUE;
31: PetscFunctionReturn(PETSC_SUCCESS);
32: }
34: static PetscErrorCode TestExplicitTS(TS ts, PetscInt order, const char subtype[])
35: {
36: PetscReal t;
37: Vec U, X, Y;
38: TSType type;
39: PetscBool done;
40: const PetscScalar *xx;
41: const PetscScalar *yy;
42: const PetscReal Tf = 1;
43: const PetscReal dt = Tf / 8;
44: const PetscReal eps = 100 * PETSC_MACHINE_EPSILON;
45: TSAdapt adapt;
46: PetscInt step;
47: TSConvergedReason reason;
49: PetscFunctionBeginUser;
50: PetscCall(TSGetType(ts, &type));
51: PetscCall(TSGetSolution(ts, &U));
52: PetscCall(VecZeroEntries(U));
53: PetscCall(TSSetStepNumber(ts, 0));
54: PetscCall(TSSetTime(ts, 0));
55: PetscCall(TSSetTimeStep(ts, dt));
56: PetscCall(TSSetMaxTime(ts, Tf));
57: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
58: PetscCall(TSSolve(ts, NULL));
59: PetscCall(TSRollBack(ts));
60: PetscCall(TSSolve(ts, NULL));
61: PetscCall(TSGetTime(ts, &t));
63: PetscCall(TSGetSolution(ts, &U));
64: PetscCall(VecDuplicate(U, &X));
65: PetscCall(TSEvaluateStep(ts, order, X, NULL));
66: PetscCall(VecGetArrayRead(X, &xx));
67: for (PetscInt i = 0; i < order; i++) {
68: PetscReal error = PetscAbsReal(PetscRealPart(xx[i]) - PetscPowReal(t, i + 1));
69: PetscCheck(error <= eps, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad solution, error %g, %s '%s'", (double)error, type, subtype);
70: }
71: PetscCall(VecRestoreArrayRead(X, &xx));
72: PetscCall(VecDestroy(&X));
74: PetscCall(TSGetSolution(ts, &U));
75: PetscCall(VecDuplicate(U, &Y));
76: PetscCall(TSEvaluateStep(ts, order - 1, Y, &done));
77: PetscCall(VecGetArrayRead(Y, &yy));
78: for (PetscInt i = 0; done && i < order - 1; i++) {
79: PetscReal error = PetscAbsReal(PetscRealPart(yy[i]) - PetscPowReal(t, i + 1));
80: PetscCheck(error <= eps, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad estimator, error %g, %s '%s'", (double)error, type, subtype);
81: }
82: PetscCall(VecRestoreArrayRead(Y, &yy));
83: PetscCall(VecDestroy(&Y));
85: PetscCall(TSGetAdapt(ts, &adapt));
86: PetscCall(TSAdaptSetCheckStage(adapt, TestCheckStage));
87: PetscCall(TSSetErrorIfStepFails(ts, PETSC_FALSE));
88: PetscCall(TSSetStepNumber(ts, 0));
89: PetscCall(TSSetTime(ts, 0));
90: PetscCall(TSSetTimeStep(ts, dt));
91: PetscCall(TSSolve(ts, NULL));
92: PetscCall(TSAdaptSetCheckStage(adapt, NULL));
93: PetscCall(TSSetErrorIfStepFails(ts, PETSC_TRUE));
94: PetscCall(TSGetStepNumber(ts, &step));
95: PetscCall(TSGetConvergedReason(ts, &reason));
96: PetscCheck(step == 2, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad step number %" PetscInt_FMT ", %s '%s'", step, type, subtype);
97: PetscCheck(reason == TS_DIVERGED_STEP_REJECTED, PETSC_COMM_SELF, PETSC_ERR_PLIB, "Bad reason %s, %s '%s'", TSConvergedReasons[reason], type, subtype);
98: PetscFunctionReturn(PETSC_SUCCESS);
99: }
101: static PetscErrorCode TestTSRK(TS ts, TSRKType type)
102: {
103: PetscInt order;
104: TSAdapt adapt;
105: PetscBool rk1, rk3, rk4;
106: TSAdaptType adapttype;
107: char savetype[32];
109: PetscFunctionBeginUser;
110: PetscCall(TSRKSetType(ts, type));
111: PetscCall(TSRKGetType(ts, &type));
112: PetscCall(TSRKGetOrder(ts, &order));
114: PetscCall(TSGetAdapt(ts, &adapt));
115: PetscCall(TSAdaptGetType(adapt, &adapttype));
116: PetscCall(PetscStrncpy(savetype, adapttype, sizeof(savetype)));
117: PetscCall(PetscStrcmp(type, TSRK1FE, &rk1));
118: PetscCall(PetscStrcmp(type, TSRK3, &rk3));
119: PetscCall(PetscStrcmp(type, TSRK4, &rk4));
120: if (rk1 || rk3 || rk4) PetscCall(TSAdaptSetType(adapt, TSADAPTNONE));
122: PetscCall(TestExplicitTS(ts, order, type));
124: PetscCall(TSGetAdapt(ts, &adapt));
125: PetscCall(TSAdaptSetType(adapt, savetype));
126: PetscFunctionReturn(PETSC_SUCCESS);
127: }
129: int main(int argc, char *argv[])
130: {
131: TS ts;
132: Vec X;
133: PetscInt N = 9;
135: PetscFunctionBeginUser;
136: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
138: PetscCall(TSCreate(PETSC_COMM_SELF, &ts));
139: PetscCall(TSSetType(ts, TSRK));
140: PetscCall(TSSetRHSFunction(ts, NULL, RHSFunction, NULL));
141: PetscCall(VecCreateSeq(PETSC_COMM_SELF, N, &X));
142: PetscCall(TSSetSolution(ts, X));
143: PetscCall(VecDestroy(&X));
144: PetscCall(TSSetFromOptions(ts));
146: PetscCall(TestTSRK(ts, TSRK1FE));
147: PetscCall(TestTSRK(ts, TSRK2A));
148: PetscCall(TestTSRK(ts, TSRK3));
149: PetscCall(TestTSRK(ts, TSRK3BS));
150: PetscCall(TestTSRK(ts, TSRK4));
151: PetscCall(TestTSRK(ts, TSRK5F));
152: PetscCall(TestTSRK(ts, TSRK5DP));
153: PetscCall(TestTSRK(ts, TSRK5BS));
154: PetscCall(TestTSRK(ts, TSRK6VR));
155: PetscCall(TestTSRK(ts, TSRK7VR));
156: PetscCall(TestTSRK(ts, TSRK8VR));
158: PetscCall(TSRollBack(ts));
159: PetscCall(TSDestroy(&ts));
161: PetscCall(PetscFinalize());
162: return 0;
163: }
165: /*TEST
167: testset:
168: output_file: output/empty.out
169: test:
170: suffix: 0
171: requires: !single
172: test:
173: suffix: 1
174: args: -ts_adapt_type none
175: test:
176: suffix: 2
177: requires: !single
178: args: -ts_adapt_type basic
179: test:
180: suffix: 3
181: args: -ts_adapt_type dsp
183: TEST*/