1 static char help[] = "Solves the trivial ODE du/dt = 1, u(0) = 0. \n\n";
2
3 #include <petscts.h>
4 #include <petscpc.h>
5
6 static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
7 static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
8
9 static PetscErrorCode PreStep(TS);
10 static PetscErrorCode PostStep(TS);
11 static PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
12 static PetscErrorCode Event(TS, PetscReal, Vec, PetscReal *, void *);
13 static PetscErrorCode PostEvent(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *);
14 static PetscErrorCode TransferSetUp(TS, PetscInt, PetscReal, Vec, PetscBool *, void *);
15 static PetscErrorCode Transfer(TS, PetscInt, Vec[], Vec[], void *);
16
main(int argc,char ** argv)17 int main(int argc, char **argv)
18 {
19 TS ts;
20 PetscInt n, ntransfer[] = {2, 2};
21 const PetscInt n_end = 11;
22 PetscReal t;
23 const PetscReal t_end = 11;
24 Vec x;
25 Vec f;
26 Mat A;
27
28 PetscFunctionBeginUser;
29 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
30
31 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
32
33 PetscCall(VecCreate(PETSC_COMM_WORLD, &f));
34 PetscCall(VecSetSizes(f, 1, PETSC_DECIDE));
35 PetscCall(VecSetFromOptions(f));
36 PetscCall(VecSetUp(f));
37 PetscCall(TSSetRHSFunction(ts, f, RHSFunction, NULL));
38 PetscCall(VecDestroy(&f));
39
40 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
41 PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE));
42 PetscCall(MatSetFromOptions(A));
43 PetscCall(MatSetUp(A));
44 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
45 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
46 /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
47 PetscCall(MatShift(A, (PetscReal)1));
48 PetscCall(MatShift(A, (PetscReal)-1));
49 PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
50 PetscCall(MatDestroy(&A));
51
52 PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
53 PetscCall(VecSetSizes(x, 1, PETSC_DECIDE));
54 PetscCall(VecSetFromOptions(x));
55 PetscCall(VecSetUp(x));
56 PetscCall(TSSetSolution(ts, x));
57 PetscCall(VecDestroy(&x));
58
59 PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
60 PetscCall(TSSetPreStep(ts, PreStep));
61 PetscCall(TSSetPostStep(ts, PostStep));
62
63 {
64 TSAdapt adapt;
65 PetscCall(TSGetAdapt(ts, &adapt));
66 PetscCall(TSAdaptSetType(adapt, TSADAPTNONE));
67 }
68 {
69 PetscInt direction[3];
70 PetscBool terminate[3];
71 direction[0] = +1;
72 terminate[0] = PETSC_FALSE;
73 direction[1] = -1;
74 terminate[1] = PETSC_FALSE;
75 direction[2] = 0;
76 terminate[2] = PETSC_FALSE;
77 PetscCall(TSSetTimeStep(ts, 1));
78 PetscCall(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL));
79 }
80 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
81 PetscCall(TSSetResize(ts, PETSC_TRUE, TransferSetUp, Transfer, ntransfer));
82 PetscCall(TSSetFromOptions(ts));
83
84 /* --- First Solve --- */
85
86 PetscCall(TSSetStepNumber(ts, 0));
87 PetscCall(TSSetTimeStep(ts, 1));
88 PetscCall(TSSetTime(ts, 0));
89 PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL));
90 PetscCall(TSSetMaxSteps(ts, 3));
91
92 PetscCall(TSGetTime(ts, &t));
93 PetscCall(TSGetSolution(ts, &x));
94 PetscCall(VecSet(x, t));
95 while (t < t_end) {
96 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
97 PetscCall(TSSolve(ts, NULL));
98 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
99 PetscCall(TSGetTime(ts, &t));
100 PetscCall(TSGetStepNumber(ts, &n));
101 PetscCall(TSSetMaxSteps(ts, PetscMin(n + 3, n_end)));
102 }
103 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
104 PetscCall(TSSolve(ts, NULL));
105 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
106
107 /* --- Second Solve --- */
108
109 PetscCall(TSSetStepNumber(ts, 0));
110 PetscCall(TSSetTimeStep(ts, 1));
111 PetscCall(TSSetTime(ts, 0));
112 PetscCall(TSSetMaxTime(ts, 3));
113 PetscCall(TSSetMaxSteps(ts, PETSC_INT_MAX));
114
115 PetscCall(TSGetTime(ts, &t));
116 PetscCall(TSGetSolution(ts, &x));
117 PetscCall(VecSet(x, t));
118 while (t < t_end) {
119 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
120 PetscCall(TSSolve(ts, NULL));
121 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
122 PetscCall(TSGetTime(ts, &t));
123 PetscCall(TSSetMaxTime(ts, PetscMin(t + 3, t_end)));
124 }
125 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
126 PetscCall(TSSolve(ts, NULL));
127 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
128
129 /* --- */
130
131 PetscCall(TSDestroy(&ts));
132
133 PetscCall(PetscFinalize());
134 return 0;
135 }
136
137 /* -------------------------------------------------------------------*/
138
RHSFunction(TS ts,PetscReal t,Vec x,Vec f,PetscCtx ctx)139 PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, PetscCtx ctx)
140 {
141 PetscFunctionBeginUser;
142 PetscCall(VecSet(f, (PetscReal)1));
143 PetscFunctionReturn(PETSC_SUCCESS);
144 }
145
RHSJacobian(TS ts,PetscReal t,Vec x,Mat A,Mat B,PetscCtx ctx)146 PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, PetscCtx ctx)
147 {
148 PetscFunctionBeginUser;
149 PetscCall(MatZeroEntries(B));
150 if (B != A) PetscCall(MatZeroEntries(A));
151 PetscFunctionReturn(PETSC_SUCCESS);
152 }
153
PreStep(TS ts)154 PetscErrorCode PreStep(TS ts)
155 {
156 PetscInt n;
157 PetscReal t;
158 Vec x;
159 const PetscScalar *a;
160 PetscBool flg;
161
162 PetscFunctionBeginUser;
163 PetscCall(TSGetStepNumber(ts, &n));
164 PetscCall(TSGetTime(ts, &t));
165 PetscCall(TSGetSolution(ts, &x));
166 PetscCall(VecGetArrayRead(x, &a));
167 PetscCall(TSGetStepResize(ts, &flg));
168 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g%s\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]), flg ? " resized" : ""));
169 PetscCall(VecRestoreArrayRead(x, &a));
170 PetscFunctionReturn(PETSC_SUCCESS);
171 }
172
PostStep(TS ts)173 PetscErrorCode PostStep(TS ts)
174 {
175 PetscInt n;
176 PetscReal t;
177 Vec x;
178 const PetscScalar *a;
179
180 PetscFunctionBeginUser;
181 PetscCall(TSGetStepNumber(ts, &n));
182 PetscCall(TSGetTime(ts, &t));
183 PetscCall(TSGetSolution(ts, &x));
184 PetscCall(VecGetArrayRead(x, &a));
185 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
186 PetscCall(VecRestoreArrayRead(x, &a));
187 PetscFunctionReturn(PETSC_SUCCESS);
188 }
189
Monitor(TS ts,PetscInt n,PetscReal t,Vec x,PetscCtx ctx)190 PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, PetscCtx ctx)
191 {
192 const PetscScalar *a;
193
194 PetscFunctionBeginUser;
195 PetscCall(VecGetArrayRead(x, &a));
196 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
197 PetscCall(VecRestoreArrayRead(x, &a));
198 PetscFunctionReturn(PETSC_SUCCESS);
199 }
200
Event(TS ts,PetscReal t,Vec x,PetscReal * fvalue,PetscCtx ctx)201 PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscReal *fvalue, PetscCtx ctx)
202 {
203 PetscFunctionBeginUser;
204 fvalue[0] = t - 5;
205 fvalue[1] = 7 - t;
206 fvalue[2] = t - 9;
207 PetscFunctionReturn(PETSC_SUCCESS);
208 }
209
PostEvent(TS ts,PetscInt nevents,PetscInt event_list[],PetscReal t,Vec x,PetscBool forwardsolve,PetscCtx ctx)210 PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, PetscCtx ctx)
211 {
212 PetscInt i;
213 const PetscScalar *a;
214
215 PetscFunctionBeginUser;
216 PetscCall(TSGetStepNumber(ts, &i));
217 PetscCall(VecGetArrayRead(x, &a));
218 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0])));
219 PetscCall(VecRestoreArrayRead(x, &a));
220 PetscFunctionReturn(PETSC_SUCCESS);
221 }
222
TransferSetUp(TS ts,PetscInt n,PetscReal t,Vec x,PetscBool * flg,PetscCtx ctx)223 PetscErrorCode TransferSetUp(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, PetscCtx ctx)
224 {
225 PetscInt *nt = (PetscInt *)ctx;
226
227 PetscFunctionBeginUser;
228 if (n == 1) {
229 nt[0] = 2;
230 nt[1] = 2;
231 }
232 *flg = (PetscBool)(nt[0] && n && n % (nt[0]) == 0);
233 if (*flg) nt[0] += nt[1];
234 PetscFunctionReturn(PETSC_SUCCESS);
235 }
236
Transfer(TS ts,PetscInt nv,Vec vin[],Vec vout[],PetscCtx ctx)237 PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], PetscCtx ctx)
238 {
239 PetscInt i;
240
241 PetscFunctionBeginUser;
242 PetscCall(TSGetStepNumber(ts, &i));
243 PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " nv %" PetscInt_FMT "\n", PETSC_FUNCTION_NAME, i, nv));
244 for (i = 0; i < nv; i++) {
245 PetscCall(VecDuplicate(vin[i], &vout[i]));
246 PetscCall(VecCopy(vin[i], vout[i]));
247 }
248 PetscFunctionReturn(PETSC_SUCCESS);
249 }
250
251 /*TEST
252
253 test:
254 suffix: euler
255 diff_args: -j
256 args: -ts_type euler
257 output_file: output/ex1.out
258
259 test:
260 suffix: ssp
261 diff_args: -j
262 args: -ts_type ssp
263 output_file: output/ex1.out
264
265 test:
266 suffix: rk
267 diff_args: -j
268 args: -ts_type rk
269 output_file: output/ex1.out
270
271 test:
272 suffix: beuler
273 diff_args: -j
274 args: -ts_type beuler
275 output_file: output/ex1_theta.out
276
277 test:
278 suffix: cn
279 diff_args: -j
280 args: -ts_type cn
281 output_file: output/ex1_theta.out
282
283 test:
284 suffix: theta
285 args: -ts_type theta
286 diff_args: -j
287 output_file: output/ex1_theta.out
288
289 test:
290 suffix: bdf
291 diff_args: -j
292 args: -ts_type bdf
293 output_file: output/ex1_bdf.out
294
295 test:
296 suffix: alpha
297 diff_args: -j
298 args: -ts_type alpha
299 output_file: output/ex1_alpha.out
300
301 test:
302 suffix: rosw
303 diff_args: -j
304 args: -ts_type rosw
305 output_file: output/ex1.out
306
307 test:
308 suffix: arkimex
309 diff_args: -j
310 args: -ts_type arkimex
311 output_file: output/ex1_arkimex.out
312
313 TEST*/
314