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;
26: PetscFunctionBeginUser;
27: PetscCall(PetscInitialize(&argc, &argv, (char *)0, help));
29: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
31: PetscCall(VecCreate(PETSC_COMM_WORLD, &f));
32: PetscCall(VecSetSizes(f, 1, PETSC_DECIDE));
33: PetscCall(VecSetFromOptions(f));
34: PetscCall(VecSetUp(f));
35: PetscCall(TSSetRHSFunction(ts, f, RHSFunction, NULL));
36: PetscCall(VecDestroy(&f));
38: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
39: PetscCall(MatSetSizes(A, 1, 1, PETSC_DECIDE, PETSC_DECIDE));
40: PetscCall(MatSetFromOptions(A));
41: PetscCall(MatSetUp(A));
42: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
43: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
44: /* ensure that the Jacobian matrix has diagonal entries since that is required by TS */
45: PetscCall(MatShift(A, (PetscReal)1));
46: PetscCall(MatShift(A, (PetscReal)-1));
47: PetscCall(TSSetRHSJacobian(ts, A, A, RHSJacobian, NULL));
48: PetscCall(MatDestroy(&A));
50: PetscCall(VecCreate(PETSC_COMM_WORLD, &x));
51: PetscCall(VecSetSizes(x, 1, PETSC_DECIDE));
52: PetscCall(VecSetFromOptions(x));
53: PetscCall(VecSetUp(x));
54: PetscCall(TSSetSolution(ts, x));
55: PetscCall(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: PetscCall(TSSetEventHandler(ts, 3, direction, terminate, Event, PostEvent, NULL));
67: }
68: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
69: PetscCall(TSSetMaxTime(ts, t_end));
70: PetscCall(TSSetFromOptions(ts));
72: PetscCall(TSSolve(ts, NULL));
74: PetscCall(TSDestroy(&ts));
75: PetscCall(PetscFinalize());
76: return 0;
77: }
79: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec x, Vec f, void *ctx)
80: {
81: PetscFunctionBeginUser;
82: PetscCall(VecSet(f, (PetscReal)1));
83: PetscFunctionReturn(PETSC_SUCCESS);
84: }
86: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec x, Mat A, Mat B, void *ctx)
87: {
88: PetscFunctionBeginUser;
89: PetscCall(MatZeroEntries(B));
90: if (B != A) PetscCall(MatZeroEntries(A));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: PetscErrorCode Event(TS ts, PetscReal t, Vec x, PetscScalar *fvalue, void *ctx)
95: {
96: PetscFunctionBeginUser;
97: fvalue[0] = t - 1.1;
98: fvalue[1] = 1.2 - t;
99: fvalue[2] = t - 1.3;
100: PetscFunctionReturn(PETSC_SUCCESS);
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;
108: PetscFunctionBeginUser;
109: PetscCall(TSGetStepNumber(ts, &i));
110: PetscCall(VecGetArrayRead(x, &a));
111: PetscCall(PetscPrintf(PetscObjectComm((PetscObject)ts), "%-10s-> step %" PetscInt_FMT " time %g value %g\n", PETSC_FUNCTION_NAME, i, (double)t, (double)PetscRealPart(a[0])));
112: PetscCall(VecRestoreArrayRead(x, &a));
113: PetscFunctionReturn(PETSC_SUCCESS);
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*/