1 static char help[] = "Adjoint and tangent linear sensitivity analysis of the basic equation for generator stability analysis.\n";
2
3 /*F
4
5 \begin{eqnarray}
6 \frac{d \theta}{dt} = \omega_b (\omega - \omega_s)
7 \frac{2 H}{\omega_s}\frac{d \omega}{dt} & = & P_m - P_max \sin(\theta) -D(\omega - \omega_s)\\
8 \end{eqnarray}
9
10 F*/
11
12 /*
13 This code demonstrate the sensitivity analysis interface to a system of ordinary differential equations with discontinuities.
14 It computes the sensitivities of an integral cost function
15 \int c*max(0,\theta(t)-u_s)^beta dt
16 w.r.t. initial conditions and the parameter P_m.
17 Backward Euler method is used for time integration.
18 The discontinuities are detected with TSEvent.
19 */
20
21 #include <petscts.h>
22 #include "ex3.h"
23
main(int argc,char ** argv)24 int main(int argc, char **argv)
25 {
26 TS ts, quadts; /* ODE integrator */
27 Vec U; /* solution will be stored here */
28 PetscMPIInt size;
29 PetscInt n = 2;
30 AppCtx ctx;
31 PetscScalar *u;
32 PetscReal du[2] = {0.0, 0.0};
33 PetscBool ensemble = PETSC_FALSE, flg1, flg2;
34 PetscReal ftime;
35 PetscInt steps;
36 PetscScalar *x_ptr, *y_ptr, *s_ptr;
37 Vec lambda[1], q, mu[1];
38 PetscInt direction[2];
39 PetscBool terminate[2];
40 Mat qgrad;
41 Mat sp; /* Forward sensitivity matrix */
42 SAMethod sa;
43
44 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
45 Initialize program
46 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
47 PetscFunctionBeginUser;
48 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
49 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
50 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
51
52 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
53 Create necessary matrix and vectors
54 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
55 PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx.Jac));
56 PetscCall(MatSetSizes(ctx.Jac, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
57 PetscCall(MatSetType(ctx.Jac, MATDENSE));
58 PetscCall(MatSetFromOptions(ctx.Jac));
59 PetscCall(MatSetUp(ctx.Jac));
60 PetscCall(MatCreateVecs(ctx.Jac, &U, NULL));
61 PetscCall(MatCreate(PETSC_COMM_WORLD, &ctx.Jacp));
62 PetscCall(MatSetSizes(ctx.Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
63 PetscCall(MatSetFromOptions(ctx.Jacp));
64 PetscCall(MatSetUp(ctx.Jacp));
65 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &ctx.DRDP));
66 PetscCall(MatSetUp(ctx.DRDP));
67 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &ctx.DRDU));
68 PetscCall(MatSetUp(ctx.DRDU));
69
70 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
71 Set runtime options
72 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
73 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
74 {
75 ctx.beta = 2;
76 ctx.c = 10000.0;
77 ctx.u_s = 1.0;
78 ctx.omega_s = 1.0;
79 ctx.omega_b = 120.0 * PETSC_PI;
80 ctx.H = 5.0;
81 PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
82 ctx.D = 5.0;
83 PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
84 ctx.E = 1.1378;
85 ctx.V = 1.0;
86 ctx.X = 0.545;
87 ctx.Pmax = ctx.E * ctx.V / ctx.X;
88 ctx.Pmax_ini = ctx.Pmax;
89 PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
90 ctx.Pm = 1.1;
91 PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
92 ctx.tf = 0.1;
93 ctx.tcl = 0.2;
94 PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
95 PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
96 PetscCall(PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL));
97 if (ensemble) {
98 ctx.tf = -1;
99 ctx.tcl = -1;
100 }
101
102 PetscCall(VecGetArray(U, &u));
103 u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
104 u[1] = 1.0;
105 PetscCall(PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1));
106 n = 2;
107 PetscCall(PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2));
108 u[0] += du[0];
109 u[1] += du[1];
110 PetscCall(VecRestoreArray(U, &u));
111 if (flg1 || flg2) {
112 ctx.tf = -1;
113 ctx.tcl = -1;
114 }
115 sa = SA_ADJ;
116 PetscCall(PetscOptionsEnum("-sa_method", "Sensitivity analysis method (adj or tlm)", "", SAMethods, (PetscEnum)sa, (PetscEnum *)&sa, NULL));
117 }
118 PetscOptionsEnd();
119
120 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
121 Create timestepping solver context
122 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
123 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
124 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
125 PetscCall(TSSetType(ts, TSBEULER));
126 PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx));
127 PetscCall(TSSetRHSJacobian(ts, ctx.Jac, ctx.Jac, (TSRHSJacobianFn *)RHSJacobian, &ctx));
128
129 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
130 Set initial conditions
131 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
132 PetscCall(TSSetSolution(ts, U));
133
134 /* Set RHS JacobianP */
135 PetscCall(TSSetRHSJacobianP(ts, ctx.Jacp, RHSJacobianP, &ctx));
136
137 PetscCall(TSCreateQuadratureTS(ts, PETSC_FALSE, &quadts));
138 PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, &ctx));
139 PetscCall(TSSetRHSJacobian(quadts, ctx.DRDU, ctx.DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, &ctx));
140 PetscCall(TSSetRHSJacobianP(quadts, ctx.DRDP, DRDPJacobianTranspose, &ctx));
141 if (sa == SA_ADJ) {
142 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
143 Save trajectory of solution so that TSAdjointSolve() may be used
144 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
145 PetscCall(TSSetSaveTrajectory(ts));
146 PetscCall(MatCreateVecs(ctx.Jac, &lambda[0], NULL));
147 PetscCall(MatCreateVecs(ctx.Jacp, &mu[0], NULL));
148 PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
149 }
150
151 if (sa == SA_TLM) {
152 PetscScalar val[2];
153 PetscInt row[] = {0, 1}, col[] = {0};
154
155 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &qgrad));
156 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 2, 1, NULL, &sp));
157 PetscCall(TSForwardSetSensitivities(ts, 1, sp));
158 PetscCall(TSForwardSetSensitivities(quadts, 1, qgrad));
159 val[0] = 1. / PetscSqrtScalar(1. - (ctx.Pm / ctx.Pmax) * (ctx.Pm / ctx.Pmax)) / ctx.Pmax;
160 val[1] = 0.0;
161 PetscCall(MatSetValues(sp, 2, row, 1, col, val, INSERT_VALUES));
162 PetscCall(MatAssemblyBegin(sp, MAT_FINAL_ASSEMBLY));
163 PetscCall(MatAssemblyEnd(sp, MAT_FINAL_ASSEMBLY));
164 }
165
166 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
167 Set solver options
168 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
169 PetscCall(TSSetMaxTime(ts, 1.0));
170 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
171 PetscCall(TSSetTimeStep(ts, 0.03125));
172 PetscCall(TSSetFromOptions(ts));
173
174 direction[0] = direction[1] = 1;
175 terminate[0] = terminate[1] = PETSC_FALSE;
176
177 PetscCall(TSSetEventHandler(ts, 2, direction, terminate, EventFunction, PostEventFunction, (void *)&ctx));
178
179 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
180 Solve nonlinear system
181 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
182 if (ensemble) {
183 for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
184 PetscCall(VecGetArray(U, &u));
185 u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
186 u[1] = ctx.omega_s;
187 u[0] += du[0];
188 u[1] += du[1];
189 PetscCall(VecRestoreArray(U, &u));
190 PetscCall(TSSetTimeStep(ts, 0.03125));
191 PetscCall(TSSolve(ts, U));
192 }
193 } else {
194 PetscCall(TSSolve(ts, U));
195 }
196 PetscCall(TSGetSolveTime(ts, &ftime));
197 PetscCall(TSGetStepNumber(ts, &steps));
198
199 if (sa == SA_ADJ) {
200 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
201 Adjoint model starts here
202 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
203 /* Set initial conditions for the adjoint integration */
204 PetscCall(VecGetArray(lambda[0], &y_ptr));
205 y_ptr[0] = 0.0;
206 y_ptr[1] = 0.0;
207 PetscCall(VecRestoreArray(lambda[0], &y_ptr));
208
209 PetscCall(VecGetArray(mu[0], &x_ptr));
210 x_ptr[0] = 0.0;
211 PetscCall(VecRestoreArray(mu[0], &x_ptr));
212
213 PetscCall(TSAdjointSolve(ts));
214
215 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n lambda: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n"));
216 PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
217 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n mu: d[Psi(tf)]/d[pm]\n"));
218 PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
219 PetscCall(TSGetCostIntegral(ts, &q));
220 PetscCall(VecGetArray(q, &x_ptr));
221 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm)));
222 PetscCall(VecRestoreArray(q, &x_ptr));
223 PetscCall(ComputeSensiP(lambda[0], mu[0], &ctx));
224 PetscCall(VecGetArray(mu[0], &x_ptr));
225 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)x_ptr[0]));
226 PetscCall(VecRestoreArray(mu[0], &x_ptr));
227 PetscCall(VecDestroy(&lambda[0]));
228 PetscCall(VecDestroy(&mu[0]));
229 }
230 if (sa == SA_TLM) {
231 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n trajectory sensitivity: d[phi(tf)]/d[pm] d[omega(tf)]/d[pm]\n"));
232 PetscCall(MatView(sp, PETSC_VIEWER_STDOUT_WORLD));
233 PetscCall(TSGetCostIntegral(ts, &q));
234 PetscCall(VecGetArray(q, &s_ptr));
235 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(s_ptr[0] - ctx.Pm)));
236 PetscCall(VecRestoreArray(q, &s_ptr));
237 PetscCall(MatDenseGetArray(qgrad, &s_ptr));
238 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n gradient=%g\n", (double)s_ptr[0]));
239 PetscCall(MatDenseRestoreArray(qgrad, &s_ptr));
240 PetscCall(MatDestroy(&qgrad));
241 PetscCall(MatDestroy(&sp));
242 }
243 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
244 Free work space. All PETSc objects should be destroyed when they are no longer needed.
245 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
246 PetscCall(MatDestroy(&ctx.Jac));
247 PetscCall(MatDestroy(&ctx.Jacp));
248 PetscCall(MatDestroy(&ctx.DRDU));
249 PetscCall(MatDestroy(&ctx.DRDP));
250 PetscCall(VecDestroy(&U));
251 PetscCall(TSDestroy(&ts));
252 PetscCall(PetscFinalize());
253 return 0;
254 }
255
256 /*TEST
257
258 build:
259 requires: !complex !single
260
261 test:
262 args: -sa_method adj -viewer_binary_skip_info -ts_type cn -pc_type lu
263
264 test:
265 suffix: 2
266 args: -sa_method tlm -ts_type cn -pc_type lu
267
268 test:
269 suffix: 3
270 args: -sa_method adj -ts_type rk -ts_rk_type 2a -ts_adapt_type dsp
271
272 TEST*/
273