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, NULL, 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_INT_MAX));
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;
160: PetscBool flg;
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: }
173: PetscErrorCode PostStep(TS ts)
174: {
175: PetscInt n;
176: PetscReal t;
177: Vec x;
178: const PetscScalar *a;
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: }
190: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
191: {
192: const PetscScalar *a;
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: }
201: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscReal *fvalue, void *ctx)
202: {
203: PetscFunctionBeginUser;
204: fvalue[0] = t - 5;
205: fvalue[1] = 7 - t;
206: fvalue[2] = t - 9;
207: PetscFunctionReturn(PETSC_SUCCESS);
208: }
210: PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx)
211: {
212: PetscInt i;
213: const PetscScalar *a;
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: }
223: PetscErrorCode TransferSetUp(TS ts, PetscInt n, PetscReal t, Vec x, PetscBool *flg, void *ctx)
224: {
225: PetscInt *nt = (PetscInt *)ctx;
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: }
237: PetscErrorCode Transfer(TS ts, PetscInt nv, Vec vin[], Vec vout[], void *ctx)
238: {
239: PetscInt i;
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: }
251: /*TEST
253: test:
254: suffix: euler
255: diff_args: -j
256: args: -ts_type euler
257: output_file: output/ex1.out
259: test:
260: suffix: ssp
261: diff_args: -j
262: args: -ts_type ssp
263: output_file: output/ex1.out
265: test:
266: suffix: rk
267: diff_args: -j
268: args: -ts_type rk
269: output_file: output/ex1.out
271: test:
272: suffix: beuler
273: diff_args: -j
274: args: -ts_type beuler
275: output_file: output/ex1_theta.out
277: test:
278: suffix: cn
279: diff_args: -j
280: args: -ts_type cn
281: output_file: output/ex1_theta.out
283: test:
284: suffix: theta
285: args: -ts_type theta
286: diff_args: -j
287: output_file: output/ex1_theta.out
289: test:
290: suffix: bdf
291: diff_args: -j
292: args: -ts_type bdf
293: output_file: output/ex1_bdf.out
295: test:
296: suffix: alpha
297: diff_args: -j
298: args: -ts_type alpha
299: output_file: output/ex1_alpha.out
301: test:
302: suffix: rosw
303: diff_args: -j
304: args: -ts_type rosw
305: output_file: output/ex1.out
307: test:
308: suffix: arkimex
309: diff_args: -j
310: args: -ts_type arkimex
311: output_file: output/ex1_arkimex.out
313: TEST*/