Actual source code: ex1.c
1: static char help[] = "Solves the trivial ODE du/dt = 1, u(0) = 0. \n\n";
3: #include <petscts.h>
4: #include <petscpc.h>
6: static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
7: static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
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, PetscScalar *, 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 *);
17: int main(int argc, char **argv)
18: {
19: TS ts;
20: PetscInt n, ntransfer = 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;
28: PetscFunctionBeginUser;
29: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
31: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
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));
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));
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));
59: PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
60: PetscCall(TSSetPreStep(ts, PreStep));
61: PetscCall(TSSetPostStep(ts, PostStep));
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(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL));
78: }
79: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
80: PetscCall(TSSetResize(ts, TransferSetUp, Transfer, &ntransfer));
81: PetscCall(TSSetFromOptions(ts));
83: /* --- First Solve --- */
85: PetscCall(TSSetStepNumber(ts, 0));
86: PetscCall(TSSetTimeStep(ts, 1));
87: PetscCall(TSSetTime(ts, 0));
88: PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL));
89: PetscCall(TSSetMaxSteps(ts, 3));
91: PetscCall(TSGetTime(ts, &t));
92: PetscCall(TSGetSolution(ts, &x));
93: PetscCall(VecSet(x, t));
94: while (t < t_end) {
95: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
96: PetscCall(TSSolve(ts, NULL));
97: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
98: PetscCall(TSGetTime(ts, &t));
99: PetscCall(TSGetStepNumber(ts, &n));
100: PetscCall(TSSetMaxSteps(ts, PetscMin(n + 3, n_end)));
101: }
102: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
103: PetscCall(TSSolve(ts, NULL));
104: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
106: /* --- Second Solve --- */
108: PetscCall(TSSetStepNumber(ts, 0));
109: PetscCall(TSSetTimeStep(ts, 1));
110: PetscCall(TSSetTime(ts, 0));
111: PetscCall(TSSetMaxTime(ts, 3));
112: PetscCall(TSSetMaxSteps(ts, PETSC_MAX_INT));
114: PetscCall(TSGetTime(ts, &t));
115: PetscCall(TSGetSolution(ts, &x));
116: PetscCall(VecSet(x, t));
117: while (t < t_end) {
118: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
119: PetscCall(TSSolve(ts, NULL));
120: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
121: PetscCall(TSGetTime(ts, &t));
122: PetscCall(TSSetMaxTime(ts, PetscMin(t + 3, t_end)));
123: }
124: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
125: PetscCall(TSSolve(ts, NULL));
126: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
128: /* --- */
130: PetscCall(TSDestroy(&ts));
132: PetscCall(PetscFinalize());
133: return 0;
134: }
136: /* -------------------------------------------------------------------*/
138: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
139: {
140: PetscFunctionBeginUser;
141: PetscCall(VecSet(f, (PetscReal)1));
142: PetscFunctionReturn(PETSC_SUCCESS);
143: }
145: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
146: {
147: PetscFunctionBeginUser;
148: PetscCall(MatZeroEntries(B));
149: if (B != A) PetscCall(MatZeroEntries(A));
150: PetscFunctionReturn(PETSC_SUCCESS);
151: }
153: PetscErrorCode PreStep(TS ts)
154: {
155: PetscInt n;
156: PetscReal t;
157: Vec x;
158: const PetscScalar *a;
160: PetscFunctionBeginUser;
161: PetscCall(TSGetStepNumber(ts, &n));
162: PetscCall(TSGetTime(ts, &t));
163: PetscCall(TSGetSolution(ts, &x));
164: PetscCall(VecGetArrayRead(x, &a));
165: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
166: PetscCall(VecRestoreArrayRead(x, &a));
167: PetscFunctionReturn(PETSC_SUCCESS);
168: }
170: PetscErrorCode PostStep(TS ts)
171: {
172: PetscInt n;
173: PetscReal t;
174: Vec x;
175: const PetscScalar *a;
177: PetscFunctionBeginUser;
178: PetscCall(TSGetStepNumber(ts, &n));
179: PetscCall(TSGetTime(ts, &t));
180: PetscCall(TSGetSolution(ts, &x));
181: PetscCall(VecGetArrayRead(x, &a));
182: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
183: PetscCall(VecRestoreArrayRead(x, &a));
184: PetscFunctionReturn(PETSC_SUCCESS);
185: }
187: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
188: {
189: const PetscScalar *a;
191: PetscFunctionBeginUser;
192: PetscCall(VecGetArrayRead(x, &a));
193: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
194: PetscCall(VecRestoreArrayRead(x, &a));
195: PetscFunctionReturn(PETSC_SUCCESS);
196: }
198: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
199: {
200: PetscFunctionBeginUser;
201: fvalue[0] = t - 5;
202: fvalue[1] = 7 - t;
203: fvalue[2] = t - 9;
204: PetscFunctionReturn(PETSC_SUCCESS);
205: }
207: PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx)
208: {
209: PetscInt i;
210: const PetscScalar *a;
212: PetscFunctionBeginUser;
213: PetscCall(TSGetStepNumber(ts, &i));
214: PetscCall(VecGetArrayRead(x, &a));
215: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0])));
216: PetscCall(VecRestoreArrayRead(x, &a));
217: PetscFunctionReturn(PETSC_SUCCESS);
218: }
220: PetscErrorCode TransferSetUp(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *ctx)
221: {
222: PetscInt *nt = (PetscInt *)ctx;
224: PetscFunctionBeginUser;
225: *flg = (PetscBool)(*nt && n && n % (*nt) == 0);
226: PetscFunctionReturn(PETSC_SUCCESS);
227: }
229: PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *ctx)
230: {
231: PetscInt i;
233: PetscFunctionBeginUser;
234: PetscCall(TSGetStepNumber(ts, &i));
235: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " nv %" PetscInt_FMT "\n", PETSC_FUNCTION_NAME, i, nv));
236: for (i = 0; i < nv; i++) {
237: PetscCall(VecDuplicate(vin[i], &vout[i]));
238: PetscCall(VecCopy(vin[i], vout[i]));
239: }
240: PetscFunctionReturn(PETSC_SUCCESS);
241: }
243: /*TEST
245: test:
246: suffix: euler
247: diff_args: -j
248: args: -ts_type euler
249: output_file: output/ex1.out
251: test:
252: suffix: ssp
253: diff_args: -j
254: args: -ts_type ssp
255: output_file: output/ex1.out
257: test:
258: suffix: rk
259: diff_args: -j
260: args: -ts_type rk
261: output_file: output/ex1.out
263: test:
264: suffix: beuler
265: diff_args: -j
266: args: -ts_type beuler
267: output_file: output/ex1.out
269: test:
270: suffix: cn
271: diff_args: -j
272: args: -ts_type cn
273: output_file: output/ex1.out
275: test:
276: suffix: theta
277: args: -ts_type theta
278: diff_args: -j
279: output_file: output/ex1.out
281: test:
282: suffix: bdf
283: diff_args: -j
284: args: -ts_type bdf
285: output_file: output/ex1_bdf.out
287: test:
288: suffix: alpha
289: diff_args: -j
290: args: -ts_type alpha
291: output_file: output/ex1.out
293: test:
294: suffix: rosw
295: diff_args: -j
296: args: -ts_type rosw
297: output_file: output/ex1.out
299: test:
300: suffix: arkimex
301: diff_args: -j
302: args: -ts_type arkimex
303: output_file: output/ex1.out
305: TEST*/