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