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, 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 *);
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;
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(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));
84: /* --- First Solve --- */
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));
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"));
107: /* --- Second Solve --- */
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_MAX_INT));
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"));
129: /* --- */
131: PetscCall(TSDestroy(&ts));
133: PetscCall(PetscFinalize());
134: return 0;
135: }
137: /* -------------------------------------------------------------------*/
139: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
140: {
141: PetscFunctionBeginUser;
142: PetscCall(VecSet(f, (PetscReal)1));
143: PetscFunctionReturn(PETSC_SUCCESS);
144: }
146: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
147: {
148: PetscFunctionBeginUser;
149: PetscCall(MatZeroEntries(B));
150: if (B != A) PetscCall(MatZeroEntries(A));
151: PetscFunctionReturn(PETSC_SUCCESS);
152: }
154: PetscErrorCode PreStep(TS ts)
155: {
156: PetscInt n;
157: PetscReal t;
158: Vec x;
159: const PetscScalar *a;
161: PetscFunctionBeginUser;
162: PetscCall(TSGetStepNumber(ts, &n));
163: PetscCall(TSGetTime(ts, &t));
164: PetscCall(TSGetSolution(ts, &x));
165: PetscCall(VecGetArrayRead(x, &a));
166: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
167: PetscCall(VecRestoreArrayRead(x, &a));
168: PetscFunctionReturn(PETSC_SUCCESS);
169: }
171: PetscErrorCode PostStep(TS ts)
172: {
173: PetscInt n;
174: PetscReal t;
175: Vec x;
176: const PetscScalar *a;
178: PetscFunctionBeginUser;
179: PetscCall(TSGetStepNumber(ts, &n));
180: PetscCall(TSGetTime(ts, &t));
181: PetscCall(TSGetSolution(ts, &x));
182: PetscCall(VecGetArrayRead(x, &a));
183: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
184: PetscCall(VecRestoreArrayRead(x, &a));
185: PetscFunctionReturn(PETSC_SUCCESS);
186: }
188: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
189: {
190: const PetscScalar *a;
192: PetscFunctionBeginUser;
193: PetscCall(VecGetArrayRead(x, &a));
194: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0])));
195: PetscCall(VecRestoreArrayRead(x, &a));
196: PetscFunctionReturn(PETSC_SUCCESS);
197: }
199: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscReal *fvalue, void *ctx)
200: {
201: PetscFunctionBeginUser;
202: fvalue[0] = t - 5;
203: fvalue[1] = 7 - t;
204: fvalue[2] = t - 9;
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx)
209: {
210: PetscInt i;
211: const PetscScalar *a;
213: PetscFunctionBeginUser;
214: PetscCall(TSGetStepNumber(ts, &i));
215: PetscCall(VecGetArrayRead(x, &a));
216: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0])));
217: PetscCall(VecRestoreArrayRead(x, &a));
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: PetscErrorCode TransferSetUp(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *ctx)
222: {
223: PetscInt *nt = (PetscInt *)ctx;
225: PetscFunctionBeginUser;
226: if (n == 1) {
227: nt[0] = 2;
228: nt[1] = 2;
229: }
230: *flg = (PetscBool)(nt[0] && n && n % (nt[0]) == 0);
231: if (*flg) nt[0] += nt[1];
232: PetscFunctionReturn(PETSC_SUCCESS);
233: }
235: PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *ctx)
236: {
237: PetscInt i;
239: PetscFunctionBeginUser;
240: PetscCall(TSGetStepNumber(ts, &i));
241: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " nv %" PetscInt_FMT "\n", PETSC_FUNCTION_NAME, i, nv));
242: for (i = 0; i < nv; i++) {
243: PetscCall(VecDuplicate(vin[i], &vout[i]));
244: PetscCall(VecCopy(vin[i], vout[i]));
245: }
246: PetscFunctionReturn(PETSC_SUCCESS);
247: }
249: /*TEST
251: test:
252: suffix: euler
253: diff_args: -j
254: args: -ts_type euler
255: output_file: output/ex1.out
257: test:
258: suffix: ssp
259: diff_args: -j
260: args: -ts_type ssp
261: output_file: output/ex1.out
263: test:
264: suffix: rk
265: diff_args: -j
266: args: -ts_type rk
267: output_file: output/ex1.out
269: test:
270: suffix: beuler
271: diff_args: -j
272: args: -ts_type beuler
273: output_file: output/ex1_theta.out
275: test:
276: suffix: cn
277: diff_args: -j
278: args: -ts_type cn
279: output_file: output/ex1_theta.out
281: test:
282: suffix: theta
283: args: -ts_type theta
284: diff_args: -j
285: output_file: output/ex1_theta.out
287: test:
288: suffix: bdf
289: diff_args: -j
290: args: -ts_type bdf
291: output_file: output/ex1_bdf.out
293: test:
294: suffix: alpha
295: diff_args: -j
296: args: -ts_type alpha
297: output_file: output/ex1_alpha.out
299: test:
300: suffix: rosw
301: diff_args: -j
302: args: -ts_type rosw
303: output_file: output/ex1.out
305: test:
306: suffix: arkimex
307: diff_args: -j
308: args: -ts_type arkimex
309: output_file: output/ex1_arkimex.out
311: TEST*/