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;
28: PetscInitialize(&argc, &argv, (char *)0, help);
30: TSCreate(PETSC_COMM_WORLD, &ts);
32: VecCreate(PETSC_COMM_WORLD, &f);
33: VecSetSizes(f, 1, PETSC_DECIDE);
34: VecSetFromOptions(f);
35: VecSetUp(f);
36: TSSetRHSFunction(ts, f, RHSFunction, NULL);
37: VecDestroy(&f);
39: MatCreate(PETSC_COMM_WORLD, &A);
40: MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE);
41: MatSetFromOptions(A);
42: MatSetUp(A);
43: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
44: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
45: /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
46: MatShift(A, (PetscReal)1);
47: MatShift(A, (PetscReal)-1);
48: TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL);
49: MatDestroy(&A);
51: VecCreate(PETSC_COMM_WORLD, &x);
52: VecSetSizes(x, 1, PETSC_DECIDE);
53: VecSetFromOptions(x);
54: VecSetUp(x);
55: TSSetSolution(ts, x);
56: VecDestroy(&x);
58: TSMonitorSet(ts, Monitor, NULL, NULL);
59: TSSetPreStep(ts, PreStep);
60: TSSetPostStep(ts, PostStep);
62: {
63: TSAdapt adapt;
64: TSGetAdapt(ts, &adapt);
65: 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: TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL);
77: }
78: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
79: TSSetFromOptions(ts);
81: /* --- First Solve --- */
83: TSSetStepNumber(ts, 0);
84: TSSetTimeStep(ts, 1);
85: TSSetTime(ts, 0);
86: TSSetMaxTime(ts, PETSC_MAX_REAL);
87: TSSetMaxSteps(ts, 3);
89: TSGetTime(ts, &t);
90: TSGetSolution(ts, &x);
91: VecSet(x, t);
92: while (t < t_end) {
93: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n");
94: TSSolve(ts, NULL);
95: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n");
96: TSGetTime(ts, &t);
97: TSGetStepNumber(ts, &n);
98: TSSetMaxSteps(ts, PetscMin(n + 3, n_end));
99: }
100: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n");
101: TSSolve(ts, NULL);
102: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n");
104: /* --- Second Solve --- */
106: TSSetStepNumber(ts, 0);
107: TSSetTimeStep(ts, 1);
108: TSSetTime(ts, 0);
109: TSSetMaxTime(ts, 3);
110: TSSetMaxSteps(ts, PETSC_MAX_INT);
112: TSGetTime(ts, &t);
113: TSGetSolution(ts, &x);
114: VecSet(x, t);
115: while (t < t_end) {
116: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n");
117: TSSolve(ts, NULL);
118: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n");
119: TSGetTime(ts, &t);
120: TSSetMaxTime(ts, PetscMin(t + 3, t_end));
121: }
122: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: Begin\n");
123: TSSolve(ts, NULL);
124: PetscPrintf(PetscObjectComm((PetscObject)ts), "TSSolve: End\n\n");
126: /* --- */
128: TSDestroy(&ts);
130: PetscFinalize();
131: return 0;
132: }
134: /* -------------------------------------------------------------------*/
136: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
137: {
139: VecSet(f, (PetscReal)1);
140: return 0;
141: }
143: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
144: {
146: MatZeroEntries(B);
147: if (B != A) MatZeroEntries(A);
148: return 0;
149: }
151: PetscErrorCode PreStep(TS ts)
152: {
153: PetscInt n;
154: PetscReal t;
155: Vec x;
156: const PetscScalar *a;
159: TSGetStepNumber(ts, &n);
160: TSGetTime(ts, &t);
161: TSGetSolution(ts, &x);
162: VecGetArrayRead(x, &a);
163: PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]));
164: VecRestoreArrayRead(x, &a);
165: return 0;
166: }
168: PetscErrorCode PostStep(TS ts)
169: {
170: PetscInt n;
171: PetscReal t;
172: Vec x;
173: const PetscScalar *a;
176: TSGetStepNumber(ts, &n);
177: TSGetTime(ts, &t);
178: TSGetSolution(ts, &x);
179: VecGetArrayRead(x, &a);
180: PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]));
181: VecRestoreArrayRead(x, &a);
182: return 0;
183: }
185: PetscErrorCode Monitor(TS ts, PetscInt n, PetscReal t, Vec x, void *ctx)
186: {
187: const PetscScalar *a;
190: VecGetArrayRead(x, &a);
191: PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, n, (double)t, (double)PetscRealPart(a[0]));
192: VecRestoreArrayRead(x, &a);
193: return 0;
194: }
196: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
197: {
199: fvalue[0] = t - 5;
200: fvalue[1] = 7 - t;
201: fvalue[2] = t - 9;
202: return 0;
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;
211: TSGetStepNumber(ts, &i);
212: VecGetArrayRead(x, &a);
213: PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0]));
214: VecRestoreArrayRead(x, &a);
215: return 0;
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*/