Actual source code: ex1fwd.c
1: static char help[] = "Trajectory sensitivity of a hybrid system with state-dependent switchings.\n";
3: /*
4: The dynamics is described by the ODE
5: u_t = A_i u
7: where A_1 = [ 1 -100
8: 10 1 ],
9: A_2 = [ 1 10
10: -100 1 ].
11: The index i changes from 1 to 2 when u[1]=2.75u[0] and from 2 to 1 when u[1]=0.36u[0].
12: Initially u=[0 1]^T and i=1.
14: References:
15: + * - H. Zhang, S. Abhyankar, E. Constantinescu, M. Mihai, Discrete Adjoint Sensitivity Analysis of Hybrid Dynamical Systems With Switching,
16: IEEE Transactions on Circuits and Systems I: Regular Papers, 64(5), May 2017
17: - * - I. A. Hiskens, M.A. Pai, Trajectory Sensitivity Analysis of Hybrid Systems, IEEE Transactions on Circuits and Systems, Vol 47, No 2, February 2000
18: */
20: #include <petscts.h>
22: typedef struct {
23: PetscScalar lambda1;
24: PetscScalar lambda2;
25: PetscInt mode; /* mode flag*/
26: PetscReal print_time;
27: } AppCtx;
29: PetscErrorCode MyMonitor(TS ts, PetscInt stepnum, PetscReal time, Vec U, void *ctx)
30: {
31: AppCtx *actx = (AppCtx *)ctx;
32: Mat sp;
33: PetscScalar *u;
34: PetscInt nump;
35: FILE *f;
37: PetscFunctionBegin;
38: if (time >= actx->print_time) {
39: actx->print_time += 1. / 256.;
40: PetscCall(TSForwardGetSensitivities(ts, &nump, &sp));
41: PetscCall(MatDenseGetColumn(sp, 2, &u));
42: f = fopen("fwd_sp.out", "a");
43: PetscCall(PetscFPrintf(PETSC_COMM_WORLD, f, "%20.15lf %20.15lf %20.15lf\n", (double)time, (double)PetscRealPart(u[0]), (double)PetscRealPart(u[1])));
44: PetscCall(MatDenseRestoreColumn(sp, &u));
45: fclose(f);
46: }
47: PetscFunctionReturn(PETSC_SUCCESS);
48: }
50: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec U, PetscReal *fvalue, void *ctx)
51: {
52: AppCtx *actx = (AppCtx *)ctx;
53: const PetscScalar *u;
55: PetscFunctionBegin;
56: PetscCall(VecGetArrayRead(U, &u));
57: if (actx->mode == 1) {
58: fvalue[0] = PetscRealPart(u[1] - actx->lambda1 * u[0]);
59: } else if (actx->mode == 2) {
60: fvalue[0] = PetscRealPart(u[1] - actx->lambda2 * u[0]);
61: }
62: PetscCall(VecRestoreArrayRead(U, &u));
63: PetscFunctionReturn(PETSC_SUCCESS);
64: }
66: PetscErrorCode ShiftGradients(TS ts, Vec U, AppCtx *actx)
67: {
68: Mat sp;
69: PetscScalar *x;
70: PetscScalar *u;
71: PetscScalar tmp[2], A1[2][2], A2[2], denorm;
72: PetscInt nump;
74: PetscFunctionBegin;
75: PetscCall(TSForwardGetSensitivities(ts, &nump, &sp));
76: PetscCall(VecGetArray(U, &u));
78: if (actx->mode == 1) {
79: denorm = -actx->lambda1 * (u[0] - 100. * u[1]) + 1. * (10. * u[0] + u[1]);
80: A1[0][0] = 110. * u[1] * (-actx->lambda1) / denorm + 1.;
81: A1[1][0] = -110. * u[0] * (-actx->lambda1) / denorm;
82: A1[0][1] = 110. * u[1] * 1. / denorm;
83: A1[1][1] = -110. * u[0] * 1. / denorm + 1.;
85: A2[0] = 110. * u[1] * (-u[0]) / denorm;
86: A2[1] = -110. * u[0] * (-u[0]) / denorm;
87: } else {
88: denorm = -actx->lambda2 * (u[0] + 10. * u[1]) + 1. * (-100. * u[0] + u[1]);
89: A1[0][0] = 110. * u[1] * (actx->lambda2) / denorm + 1;
90: A1[1][0] = -110. * u[0] * (actx->lambda2) / denorm;
91: A1[0][1] = -110. * u[1] * 1. / denorm;
92: A1[1][1] = 110. * u[0] * 1. / denorm + 1.;
94: A2[0] = 0;
95: A2[1] = 0;
96: }
98: PetscCall(VecRestoreArray(U, &u));
100: PetscCall(MatDenseGetColumn(sp, 0, &x));
101: tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
102: tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
103: x[0] = tmp[0];
104: x[1] = tmp[1];
105: PetscCall(MatDenseRestoreColumn(sp, &x));
107: PetscCall(MatDenseGetColumn(sp, 1, &x));
108: tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
109: tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
110: x[0] = tmp[0];
111: x[1] = tmp[1];
112: PetscCall(MatDenseRestoreColumn(sp, &x));
114: PetscCall(MatDenseGetColumn(sp, 2, &x));
115: tmp[0] = A1[0][0] * x[0] + A1[0][1] * x[1];
116: tmp[1] = A1[1][0] * x[0] + A1[1][1] * x[1];
117: x[0] = tmp[0] + A2[0];
118: x[1] = tmp[1] + A2[1];
119: PetscCall(MatDenseRestoreColumn(sp, &x));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec U, PetscBool forwardsolve, void *ctx)
124: {
125: AppCtx *actx = (AppCtx *)ctx;
127: PetscFunctionBegin;
128: /* PetscCall(VecView(U,PETSC_VIEWER_STDOUT_WORLD)); */
129: PetscCall(ShiftGradients(ts, U, actx));
130: if (actx->mode == 1) {
131: actx->mode = 2;
132: /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"Change from mode 1 to 2 at t = %f \n",(double)t)); */
133: } else if (actx->mode == 2) {
134: actx->mode = 1;
135: /* PetscCall(PetscPrintf(PETSC_COMM_SELF,"Change from mode 2 to 1 at t = %f \n",(double)t)); */
136: }
137: PetscFunctionReturn(PETSC_SUCCESS);
138: }
140: /*
141: Defines the ODE passed to the ODE solver
142: */
143: static PetscErrorCode IFunction(TS ts, PetscReal t, Vec U, Vec Udot, Vec F, void *ctx)
144: {
145: AppCtx *actx = (AppCtx *)ctx;
146: PetscScalar *f;
147: const PetscScalar *u, *udot;
149: PetscFunctionBegin;
150: /* The next three lines allow us to access the entries of the vectors directly */
151: PetscCall(VecGetArrayRead(U, &u));
152: PetscCall(VecGetArrayRead(Udot, &udot));
153: PetscCall(VecGetArray(F, &f));
155: if (actx->mode == 1) {
156: f[0] = udot[0] - u[0] + 100 * u[1];
157: f[1] = udot[1] - 10 * u[0] - u[1];
158: } else if (actx->mode == 2) {
159: f[0] = udot[0] - u[0] - 10 * u[1];
160: f[1] = udot[1] + 100 * u[0] - u[1];
161: }
163: PetscCall(VecRestoreArrayRead(U, &u));
164: PetscCall(VecRestoreArrayRead(Udot, &udot));
165: PetscCall(VecRestoreArray(F, &f));
166: PetscFunctionReturn(PETSC_SUCCESS);
167: }
169: /*
170: Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
171: */
172: static PetscErrorCode IJacobian(TS ts, PetscReal t, Vec U, Vec Udot, PetscReal a, Mat A, Mat B, void *ctx)
173: {
174: AppCtx *actx = (AppCtx *)ctx;
175: PetscInt rowcol[] = {0, 1};
176: PetscScalar J[2][2];
177: const PetscScalar *u, *udot;
179: PetscFunctionBegin;
180: PetscCall(VecGetArrayRead(U, &u));
181: PetscCall(VecGetArrayRead(Udot, &udot));
183: if (actx->mode == 1) {
184: J[0][0] = a - 1;
185: J[0][1] = 100;
186: J[1][0] = -10;
187: J[1][1] = a - 1;
188: } else if (actx->mode == 2) {
189: J[0][0] = a - 1;
190: J[0][1] = -10;
191: J[1][0] = 100;
192: J[1][1] = a - 1;
193: }
194: PetscCall(MatSetValues(B, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
196: PetscCall(VecRestoreArrayRead(U, &u));
197: PetscCall(VecRestoreArrayRead(Udot, &udot));
199: PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
200: PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
201: if (A != B) {
202: PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
203: PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
204: }
205: PetscFunctionReturn(PETSC_SUCCESS);
206: }
208: /* Matrix JacobianP is constant so that it only needs to be evaluated once */
209: static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat Ap, void *ctx)
210: {
211: PetscFunctionBeginUser;
212: PetscFunctionReturn(PETSC_SUCCESS);
213: }
215: int main(int argc, char **argv)
216: {
217: TS ts; /* ODE integrator */
218: Vec U; /* solution will be stored here */
219: Mat A; /* Jacobian matrix */
220: Mat Ap; /* Jacobian dfdp */
221: PetscMPIInt size;
222: PetscInt n = 2;
223: PetscScalar *u;
224: AppCtx app;
225: PetscInt direction[1];
226: PetscBool terminate[1];
227: Mat sp;
228: PetscReal tend;
229: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
230: Initialize program
231: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
232: PetscFunctionBeginUser;
233: PetscCall(PetscInitialize(&argc, &argv, NULL, help));
234: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
235: PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
236: app.mode = 1;
237: app.lambda1 = 2.75;
238: app.lambda2 = 0.36;
239: app.print_time = 1. / 256.;
240: tend = 0.125;
241: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "ex1fwd options", "");
242: {
243: PetscCall(PetscOptionsReal("-lambda1", "", "", app.lambda1, &app.lambda1, NULL));
244: PetscCall(PetscOptionsReal("-lambda2", "", "", app.lambda2, &app.lambda2, NULL));
245: PetscCall(PetscOptionsReal("-tend", "", "", tend, &tend, NULL));
246: }
247: PetscOptionsEnd();
249: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
250: Create necessary matrix and vectors
251: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
252: PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
253: PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
254: PetscCall(MatSetType(A, MATDENSE));
255: PetscCall(MatSetFromOptions(A));
256: PetscCall(MatSetUp(A));
258: PetscCall(MatCreateVecs(A, &U, NULL));
260: PetscCall(MatCreate(PETSC_COMM_WORLD, &Ap));
261: PetscCall(MatSetSizes(Ap, n, 3, PETSC_DETERMINE, PETSC_DETERMINE));
262: PetscCall(MatSetType(Ap, MATDENSE));
263: PetscCall(MatSetFromOptions(Ap));
264: PetscCall(MatSetUp(Ap));
265: PetscCall(MatZeroEntries(Ap));
267: PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, n, 3, NULL, &sp));
268: PetscCall(MatZeroEntries(sp));
269: PetscCall(MatShift(sp, 1.0));
271: PetscCall(VecGetArray(U, &u));
272: u[0] = 0;
273: u[1] = 1;
274: PetscCall(VecRestoreArray(U, &u));
276: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
277: Create timestepping solver context
278: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
279: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
280: PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
281: PetscCall(TSSetType(ts, TSCN));
282: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &app));
283: PetscCall(TSSetIJacobian(ts, A, A, (TSIJacobianFn *)IJacobian, &app));
285: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
286: Set initial conditions
287: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
288: PetscCall(TSSetSolution(ts, U));
290: PetscCall(TSForwardSetSensitivities(ts, 3, sp));
291: /* Set RHS JacobianP */
292: PetscCall(TSSetRHSJacobianP(ts, Ap, RHSJacobianP, &app));
294: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
295: Set solver options
296: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
297: PetscCall(TSSetMaxTime(ts, tend));
298: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
299: PetscCall(TSSetTimeStep(ts, 1. / 256.));
300: PetscCall(TSMonitorSet(ts, MyMonitor, &app, NULL));
301: PetscCall(TSSetFromOptions(ts));
303: /* Set directions and terminate flags for the two events */
304: direction[0] = 0;
305: terminate[0] = PETSC_FALSE;
306: PetscCall(TSSetEventHandler(ts, 1, direction, terminate, EventFunction, PostEventFunction, (void *)&app));
307: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
308: Run timestepping solver
309: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
310: PetscCall(TSSolve(ts, U));
312: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
313: Free work space. All PETSc objects should be destroyed when they are no longer needed.
314: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
315: PetscCall(MatDestroy(&A));
316: PetscCall(VecDestroy(&U));
317: PetscCall(TSDestroy(&ts));
319: PetscCall(MatDestroy(&Ap));
320: PetscCall(MatDestroy(&sp));
321: PetscCall(PetscFinalize());
322: return 0;
323: }
325: /*TEST
327: build:
328: requires: !complex
330: test:
331: args: -ts_monitor
333: TEST*/