Actual source code: ex1.c
2: static char help[] = "Solves the trivial ODE du/dt = 1, u(0) = 0. \n\n";
4: #include <petscts.h>
5: #include <petscpc.h>
7: static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
8: static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
10: static PetscErrorCode PreStep(TS);
11: static PetscErrorCode PostStep(TS);
12: static PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *);
13: static PetscErrorCode Event(TS, PetscReal, Vec, PetscScalar *, void *);
14: static PetscErrorCode PostEvent(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *);
16: int main(int argc, char **argv)
17: {
18: TS ts;
19: PetscInt n;
20: const PetscInt n_end = 11;
21: PetscReal t;
22: const PetscReal t_end = 11;
23: Vec x;
24: Vec f;
25: Mat A;
27: PetscFunctionBeginUser;
28: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
30: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
32: PetscCall(VecCreate(PETSC_COMM_WORLD, &f));
33: PetscCall(VecSetSizes(f, 1, PETSC_DECIDE));
34: PetscCall(VecSetFromOptions(f));
35: PetscCall(VecSetUp(f));
36: PetscCall(TSSetRHSFunction(ts, f, RHSFunction, NULL));
37: PetscCall(VecDestroy(&f));
39: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
40: PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE));
41: PetscCall(MatSetFromOptions(A));
42: PetscCall(MatSetUp(A));
43: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
44: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
45: /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
46: PetscCall(MatShift(A, (PetscReal)1));
47: PetscCall(MatShift(A, (PetscReal)-1));
48: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
49: PetscCall(MatDestroy(&A));
51: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
52: PetscCall(VecSetSizes(x, 1, PETSC_DECIDE));
53: PetscCall(VecSetFromOptions(x));
54: PetscCall(VecSetUp(x));
55: PetscCall(TSSetSolution(ts, x));
56: PetscCall(VecDestroy(&x));
58: PetscCall(TSMonitorSet(ts, Monitor, NULL, NULL));
59: PetscCall(TSSetPreStep(ts, PreStep));
60: PetscCall(TSSetPostStep(ts, PostStep));
62: {
63: TSAdapt adapt;
64: PetscCall(TSGetAdapt(ts, &adapt));
65: PetscCall(TSAdaptSetType(adapt, TSADAPTNONE));
66: }
67: {
68: PetscInt direction[3];
69: PetscBool terminate[3];
70: direction[0] = +1;
71: terminate[0] = PETSC_FALSE;
72: direction[1] = -1;
73: terminate[1] = PETSC_FALSE;
74: direction[2] = 0;
75: terminate[2] = PETSC_FALSE;
76: PetscCall(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL));
77: }
78: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
79: PetscCall(TSSetFromOptions(ts));
81: /* --- First Solve --- */
83: PetscCall(TSSetStepNumber(ts, 0));
84: PetscCall(TSSetTimeStep(ts, 1));
85: PetscCall(TSSetTime(ts, 0));
86: PetscCall(TSSetMaxTime(ts, PETSC_MAX_REAL));
87: PetscCall(TSSetMaxSteps(ts, 3));
89: PetscCall(TSGetTime(ts, &t));
90: PetscCall(TSGetSolution(ts, &x));
91: PetscCall(VecSet(x, t));
92: while (t < t_end) {
93: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
94: PetscCall(TSSolve(ts, NULL));
95: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
96: PetscCall(TSGetTime(ts, &t));
97: PetscCall(TSGetStepNumber(ts, &n));
98: PetscCall(TSSetMaxSteps(ts, PetscMin(n + 3, n_end)));
99: }
100: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
101: PetscCall(TSSolve(ts, NULL));
102: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
104: /* --- Second Solve --- */
106: PetscCall(TSSetStepNumber(ts, 0));
107: PetscCall(TSSetTimeStep(ts, 1));
108: PetscCall(TSSetTime(ts, 0));
109: PetscCall(TSSetMaxTime(ts, 3));
110: PetscCall(TSSetMaxSteps(ts, PETSC_MAX_INT));
112: PetscCall(TSGetTime(ts, &t));
113: PetscCall(TSGetSolution(ts, &x));
114: PetscCall(VecSet(x, t));
115: while (t < t_end) {
116: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
117: PetscCall(TSSolve(ts, NULL));
118: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
119: PetscCall(TSGetTime(ts, &t));
120: PetscCall(TSSetMaxTime(ts, PetscMin(t + 3, t_end)));
121: }
122: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n"));
123: PetscCall(TSSolve(ts, NULL));
124: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n"));
126: /* --- */
128: PetscCall(TSDestroy(&ts));
130: PetscCall(PetscFinalize());
131: return 0;
132: }
134: /* -------------------------------------------------------------------*/
136: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
137: {
138: PetscFunctionBeginUser;
139: PetscCall(VecSet(f, (PetscReal)1));
140: PetscFunctionReturn(PETSC_SUCCESS);
141: }
143: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
144: {
145: PetscFunctionBeginUser;
146: PetscCall(MatZeroEntries(B));
147: if (B != A) PetscCall(MatZeroEntries(A));
148: PetscFunctionReturn(PETSC_SUCCESS);
149: }
151: PetscErrorCode PreStep(TS ts)
152: {
153: PetscInt n;
154: PetscReal t;
155: Vec x;
156: const PetscScalar *a;
158: PetscFunctionBeginUser;
159: PetscCall(TSGetStepNumber(ts, &n));
160: PetscCall(TSGetTime(ts, &t));
161: PetscCall(TSGetSolution(ts, &x));
162: PetscCall(VecGetArrayRead(x, &a));
163: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
164: PetscCall(VecRestoreArrayRead(x, &a));
165: PetscFunctionReturn(PETSC_SUCCESS);
166: }
168: PetscErrorCode PostStep(TS ts)
169: {
170: PetscInt n;
171: PetscReal t;
172: Vec x;
173: const PetscScalar *a;
175: PetscFunctionBeginUser;
176: PetscCall(TSGetStepNumber(ts, &n));
177: PetscCall(TSGetTime(ts, &t));
178: PetscCall(TSGetSolution(ts, &x));
179: PetscCall(VecGetArrayRead(x, &a));
180: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
181: PetscCall(VecRestoreArrayRead(x, &a));
182: PetscFunctionReturn(PETSC_SUCCESS);
183: }
185: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
186: {
187: const PetscScalar *a;
189: PetscFunctionBeginUser;
190: PetscCall(VecGetArrayRead(x, &a));
191: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
192: PetscCall(VecRestoreArrayRead(x, &a));
193: PetscFunctionReturn(PETSC_SUCCESS);
194: }
196: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
197: {
198: PetscFunctionBeginUser;
199: fvalue[0] = t - 5;
200: fvalue[1] = 7 - t;
201: fvalue[2] = t - 9;
202: PetscFunctionReturn(PETSC_SUCCESS);
203: }
205: PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx)
206: {
207: PetscInt i;
208: const PetscScalar *a;
210: PetscFunctionBeginUser;
211: PetscCall(TSGetStepNumber(ts, &i));
212: PetscCall(VecGetArrayRead(x, &a));
213: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0])));
214: PetscCall(VecRestoreArrayRead(x, &a));
215: PetscFunctionReturn(PETSC_SUCCESS);
216: }
218: /*TEST
220: test:
221: suffix: euler
222: args: -ts_type euler
223: output_file: output/ex1.out
225: test:
226: suffix: ssp
227: args: -ts_type ssp
228: output_file: output/ex1.out
230: test:
231: suffix: rk
232: args: -ts_type rk
233: output_file: output/ex1.out
235: test:
236: suffix: beuler
237: args: -ts_type beuler
238: output_file: output/ex1.out
240: test:
241: suffix: cn
242: args: -ts_type cn
243: output_file: output/ex1.out
245: test:
246: suffix: theta
247: args: -ts_type theta
248: output_file: output/ex1.out
250: test:
251: suffix: bdf
252: args: -ts_type bdf
253: output_file: output/ex1.out
255: test:
256: suffix: alpha
257: args: -ts_type alpha
258: output_file: output/ex1.out
260: test:
261: suffix: rosw
262: args: -ts_type rosw
263: output_file: output/ex1.out
265: test:
266: suffix: arkimex
267: args: -ts_type arkimex
268: output_file: output/ex1.out
270: TEST*/