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