Actual source code: ex16.c
2: static char help[] = "Solves the trivial ODE du/dt = 1, u(0) = 0. \n\n";
3: /*
4: This example tests TSEvent's capability to handle complicated cases.
5: Test 1: an event close to endpoint of a time step should not be detected twice.
6: Test 2: two events in the same time step should be detected in the correct order.
7: Test 3: three events in the same time step should be detected in the correct order.
8: */
10: #include <petscts.h>
12: static PetscErrorCode RHSFunction(TS, PetscReal, Vec, Vec, void *);
13: static PetscErrorCode RHSJacobian(TS, PetscReal, Vec, Mat, Mat, void *);
15: static PetscErrorCode Event(TS, PetscReal, Vec, PetscScalar *, void *);
16: static PetscErrorCode PostEvent(TS, PetscInt, PetscInt[], PetscReal, Vec, PetscBool, void *);
18: int main(int argc, char **argv)
19: {
20: TS ts;
21: const PetscReal t_end = 2.0;
22: Vec x;
23: Vec f;
24: Mat A;
27: PetscInitialize(&argc, &argv, (char *)0, help);
29: TSCreate(PETSC_COMM_WORLD, &ts);
31: VecCreate(PETSC_COMM_WORLD, &f);
32: VecSetSizes(f, 1, PETSC_DECIDE);
33: VecSetFromOptions(f);
34: VecSetUp(f);
35: TSSetRHSFunction(ts, f, RHSFunction, NULL);
36: VecDestroy(&f);
38: MatCreate(PETSC_COMM_WORLD, &A);
39: MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE);
40: MatSetFromOptions(A);
41: MatSetUp(A);
42: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
43: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
44: /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
45: MatShift(A, (PetscReal)1);
46: MatShift(A, (PetscReal)-1);
47: TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL);
48: MatDestroy(&A);
50: VecCreate(PETSC_COMM_WORLD, &x);
51: VecSetSizes(x, 1, PETSC_DECIDE);
52: VecSetFromOptions(x);
53: VecSetUp(x);
54: TSSetSolution(ts, x);
55: VecDestroy(&x);
57: {
58: PetscInt direction[3];
59: PetscBool terminate[3];
60: direction[0] = +1;
61: terminate[0] = PETSC_FALSE;
62: direction[1] = -1;
63: terminate[1] = PETSC_FALSE;
64: direction[2] = 0;
65: terminate[2] = PETSC_FALSE;
66: TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL);
67: }
68: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
69: TSSetMaxTime(ts, t_end);
70: TSSetFromOptions(ts);
72: TSSolve(ts, NULL);
74: TSDestroy(&ts);
75: PetscFinalize();
76: return 0;
77: }
79: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
80: {
82: VecSet(f, (PetscReal)1);
83: return 0;
84: }
86: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
87: {
89: MatZeroEntries(B);
90: if (B != A) MatZeroEntries(A);
91: return 0;
92: }
94: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
95: {
97: fvalue[0] = t - 1.1;
98: fvalue[1] = 1.2 - t;
99: fvalue[2] = t - 1.3;
100: return 0;
101: }
103: PetscErrorCode PostEvent(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec x, PetscBool forwardsolve, void *ctx)
104: {
105: PetscInt i;
106: const PetscScalar *a;
109: TSGetStepNumber(ts, &i);
110: VecGetArrayRead(x, &a);
111: PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0]));
112: VecRestoreArrayRead(x, &a);
113: return 0;
114: }
116: /*TEST
118: test:
119: args: -ts_type beuler -ts_dt 0.1 -ts_event_monitor
121: test:
122: suffix: 2
123: args: -ts_type beuler -ts_dt 0.2 -ts_event_monitor
125: test:
126: suffix: 3
127: args: -ts_type beuler -ts_dt 0.5 -ts_event_monitor
128: TEST*/