1 static char help[] = "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 Ensemble of initial conditions
11 ./ex9 -ensemble -ts_monitor_draw_solution_phase -1,-3,3,3 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly
12
13 Fault at .1 seconds
14 ./ex9 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly
15
16 Initial conditions same as when fault is ended
17 ./ex9 -u 0.496792,1.00932 -ts_monitor_draw_solution_phase .42,.95,.6,1.05 -ts_adapt_dt_max .01 -ts_monitor -ts_type rk -pc_type lu -ksp_type preonly
18
19 F*/
20
21 /*
22 Include "petscts.h" so that we can use TS solvers. Note that this
23 file automatically includes:
24 petscsys.h - base PETSc routines petscvec.h - vectors
25 petscmat.h - matrices
26 petscis.h - index sets petscksp.h - Krylov subspace methods
27 petscviewer.h - viewers petscpc.h - preconditioners
28 petscksp.h - linear solvers
29 */
30
31 #include <petscts.h>
32
33 typedef struct {
34 PetscScalar H, D, omega_b, omega_s, Pmax, Pm, E, V, X, u_s, c;
35 PetscInt beta;
36 PetscReal tf, tcl;
37 } AppCtx;
38
PostStepFunction(TS ts)39 PetscErrorCode PostStepFunction(TS ts)
40 {
41 Vec U;
42 PetscReal t;
43 const PetscScalar *u;
44
45 PetscFunctionBegin;
46 PetscCall(TSGetTime(ts, &t));
47 PetscCall(TSGetSolution(ts, &U));
48 PetscCall(VecGetArrayRead(U, &u));
49 PetscCall(PetscPrintf(PETSC_COMM_SELF, "delta(%3.2f) = %8.7f\n", (double)t, (double)u[0]));
50 PetscCall(VecRestoreArrayRead(U, &u));
51 PetscFunctionReturn(PETSC_SUCCESS);
52 }
53
54 /*
55 Defines the ODE passed to the ODE solver
56 */
RHSFunction(TS ts,PetscReal t,Vec U,Vec F,AppCtx * ctx)57 static PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec U, Vec F, AppCtx *ctx)
58 {
59 PetscScalar *f, Pmax;
60 const PetscScalar *u;
61
62 PetscFunctionBegin;
63 /* The next three lines allow us to access the entries of the vectors directly */
64 PetscCall(VecGetArrayRead(U, &u));
65 PetscCall(VecGetArray(F, &f));
66 if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
67 else Pmax = ctx->Pmax;
68
69 f[0] = ctx->omega_b * (u[1] - ctx->omega_s);
70 f[1] = (-Pmax * PetscSinScalar(u[0]) - ctx->D * (u[1] - ctx->omega_s) + ctx->Pm) * ctx->omega_s / (2.0 * ctx->H);
71
72 PetscCall(VecRestoreArrayRead(U, &u));
73 PetscCall(VecRestoreArray(F, &f));
74 PetscFunctionReturn(PETSC_SUCCESS);
75 }
76
77 /*
78 Defines the Jacobian of the ODE passed to the ODE solver. See TSSetIJacobian() for the meaning of a and the Jacobian.
79 */
RHSJacobian(TS ts,PetscReal t,Vec U,Mat A,Mat B,AppCtx * ctx)80 static PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec U, Mat A, Mat B, AppCtx *ctx)
81 {
82 PetscInt rowcol[] = {0, 1};
83 PetscScalar J[2][2], Pmax;
84 const PetscScalar *u;
85
86 PetscFunctionBegin;
87 PetscCall(VecGetArrayRead(U, &u));
88 if ((t > ctx->tf) && (t < ctx->tcl)) Pmax = 0.0; /* A short-circuit on the generator terminal that drives the electrical power output (Pmax*sin(delta)) to 0 */
89 else Pmax = ctx->Pmax;
90
91 J[0][0] = 0;
92 J[0][1] = ctx->omega_b;
93 J[1][1] = -ctx->D * ctx->omega_s / (2.0 * ctx->H);
94 J[1][0] = -Pmax * PetscCosScalar(u[0]) * ctx->omega_s / (2.0 * ctx->H);
95
96 PetscCall(MatSetValues(A, 2, rowcol, 2, rowcol, &J[0][0], INSERT_VALUES));
97 PetscCall(VecRestoreArrayRead(U, &u));
98
99 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
100 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
101 if (A != B) {
102 PetscCall(MatAssemblyBegin(B, MAT_FINAL_ASSEMBLY));
103 PetscCall(MatAssemblyEnd(B, MAT_FINAL_ASSEMBLY));
104 }
105 PetscFunctionReturn(PETSC_SUCCESS);
106 }
107
RHSJacobianP(TS ts,PetscReal t,Vec X,Mat A,PetscCtx ctx0)108 static PetscErrorCode RHSJacobianP(TS ts, PetscReal t, Vec X, Mat A, PetscCtx ctx0)
109 {
110 PetscInt row[] = {0, 1}, col[] = {0};
111 PetscScalar J[2][1];
112 AppCtx *ctx = (AppCtx *)ctx0;
113
114 PetscFunctionBeginUser;
115 J[0][0] = 0;
116 J[1][0] = ctx->omega_s / (2.0 * ctx->H);
117 PetscCall(MatSetValues(A, 2, row, 1, col, &J[0][0], INSERT_VALUES));
118 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
119 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
120 PetscFunctionReturn(PETSC_SUCCESS);
121 }
122
CostIntegrand(TS ts,PetscReal t,Vec U,Vec R,AppCtx * ctx)123 static PetscErrorCode CostIntegrand(TS ts, PetscReal t, Vec U, Vec R, AppCtx *ctx)
124 {
125 PetscScalar *r;
126 const PetscScalar *u;
127
128 PetscFunctionBegin;
129 PetscCall(VecGetArrayRead(U, &u));
130 PetscCall(VecGetArray(R, &r));
131 r[0] = ctx->c * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta);
132 PetscCall(VecRestoreArray(R, &r));
133 PetscCall(VecRestoreArrayRead(U, &u));
134 PetscFunctionReturn(PETSC_SUCCESS);
135 }
136
DRDUJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDU,Mat B,AppCtx * ctx)137 static PetscErrorCode DRDUJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDU, Mat B, AppCtx *ctx)
138 {
139 PetscScalar ru[1];
140 const PetscScalar *u;
141 PetscInt row[] = {0}, col[] = {0};
142
143 PetscFunctionBegin;
144 PetscCall(VecGetArrayRead(U, &u));
145 ru[0] = ctx->c * ctx->beta * PetscPowScalarInt(PetscMax(0., u[0] - ctx->u_s), ctx->beta - 1);
146 PetscCall(VecRestoreArrayRead(U, &u));
147 PetscCall(MatSetValues(DRDU, 1, row, 1, col, ru, INSERT_VALUES));
148 PetscCall(MatAssemblyBegin(DRDU, MAT_FINAL_ASSEMBLY));
149 PetscCall(MatAssemblyEnd(DRDU, MAT_FINAL_ASSEMBLY));
150 PetscFunctionReturn(PETSC_SUCCESS);
151 }
152
DRDPJacobianTranspose(TS ts,PetscReal t,Vec U,Mat DRDP,AppCtx * ctx)153 static PetscErrorCode DRDPJacobianTranspose(TS ts, PetscReal t, Vec U, Mat DRDP, AppCtx *ctx)
154 {
155 PetscFunctionBegin;
156 PetscCall(MatZeroEntries(DRDP));
157 PetscCall(MatAssemblyBegin(DRDP, MAT_FINAL_ASSEMBLY));
158 PetscCall(MatAssemblyEnd(DRDP, MAT_FINAL_ASSEMBLY));
159 PetscFunctionReturn(PETSC_SUCCESS);
160 }
161
ComputeSensiP(Vec lambda,Vec mu,AppCtx * ctx)162 PetscErrorCode ComputeSensiP(Vec lambda, Vec mu, AppCtx *ctx)
163 {
164 PetscScalar sensip;
165 const PetscScalar *x, *y;
166
167 PetscFunctionBegin;
168 PetscCall(VecGetArrayRead(lambda, &x));
169 PetscCall(VecGetArrayRead(mu, &y));
170 sensip = 1. / PetscSqrtScalar(1. - (ctx->Pm / ctx->Pmax) * (ctx->Pm / ctx->Pmax)) / ctx->Pmax * x[0] + y[0];
171 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt parameter pm: %.7f \n", (double)sensip));
172 PetscCall(VecRestoreArrayRead(lambda, &x));
173 PetscCall(VecRestoreArrayRead(mu, &y));
174 PetscFunctionReturn(PETSC_SUCCESS);
175 }
176
main(int argc,char ** argv)177 int main(int argc, char **argv)
178 {
179 TS ts, quadts; /* ODE integrator */
180 Vec U; /* solution will be stored here */
181 Mat A; /* Jacobian matrix */
182 Mat Jacp; /* Jacobian matrix */
183 Mat DRDU, DRDP;
184 PetscMPIInt size;
185 PetscInt n = 2;
186 AppCtx ctx;
187 PetscScalar *u;
188 PetscReal du[2] = {0.0, 0.0};
189 PetscBool ensemble = PETSC_FALSE, flg1, flg2;
190 PetscReal ftime;
191 PetscInt steps;
192 PetscScalar *x_ptr, *y_ptr;
193 Vec lambda[1], q, mu[1];
194
195 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
196 Initialize program
197 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
198 PetscFunctionBeginUser;
199 PetscCall(PetscInitialize(&argc, &argv, NULL, help));
200 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
201 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
202
203 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
204 Create necessary matrix and vectors
205 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
206 PetscCall(MatCreate(PETSC_COMM_WORLD, &A));
207 PetscCall(MatSetSizes(A, n, n, PETSC_DETERMINE, PETSC_DETERMINE));
208 PetscCall(MatSetType(A, MATDENSE));
209 PetscCall(MatSetFromOptions(A));
210 PetscCall(MatSetUp(A));
211
212 PetscCall(MatCreateVecs(A, &U, NULL));
213
214 PetscCall(MatCreate(PETSC_COMM_WORLD, &Jacp));
215 PetscCall(MatSetSizes(Jacp, PETSC_DECIDE, PETSC_DECIDE, 2, 1));
216 PetscCall(MatSetFromOptions(Jacp));
217 PetscCall(MatSetUp(Jacp));
218
219 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 1, NULL, &DRDP));
220 PetscCall(MatSetUp(DRDP));
221 PetscCall(MatCreateDense(PETSC_COMM_WORLD, PETSC_DECIDE, PETSC_DECIDE, 1, 2, NULL, &DRDU));
222 PetscCall(MatSetUp(DRDU));
223
224 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
225 Set runtime options
226 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
227 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Swing equation options", "");
228 {
229 ctx.beta = 2;
230 ctx.c = 10000.0;
231 ctx.u_s = 1.0;
232 ctx.omega_s = 1.0;
233 ctx.omega_b = 120.0 * PETSC_PI;
234 ctx.H = 5.0;
235 PetscCall(PetscOptionsScalar("-Inertia", "", "", ctx.H, &ctx.H, NULL));
236 ctx.D = 5.0;
237 PetscCall(PetscOptionsScalar("-D", "", "", ctx.D, &ctx.D, NULL));
238 ctx.E = 1.1378;
239 ctx.V = 1.0;
240 ctx.X = 0.545;
241 ctx.Pmax = ctx.E * ctx.V / ctx.X;
242 PetscCall(PetscOptionsScalar("-Pmax", "", "", ctx.Pmax, &ctx.Pmax, NULL));
243 ctx.Pm = 1.1;
244 PetscCall(PetscOptionsScalar("-Pm", "", "", ctx.Pm, &ctx.Pm, NULL));
245 ctx.tf = 0.1;
246 ctx.tcl = 0.2;
247 PetscCall(PetscOptionsReal("-tf", "Time to start fault", "", ctx.tf, &ctx.tf, NULL));
248 PetscCall(PetscOptionsReal("-tcl", "Time to end fault", "", ctx.tcl, &ctx.tcl, NULL));
249 PetscCall(PetscOptionsBool("-ensemble", "Run ensemble of different initial conditions", "", ensemble, &ensemble, NULL));
250 if (ensemble) {
251 ctx.tf = -1;
252 ctx.tcl = -1;
253 }
254
255 PetscCall(VecGetArray(U, &u));
256 u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
257 u[1] = 1.0;
258 PetscCall(PetscOptionsRealArray("-u", "Initial solution", "", u, &n, &flg1));
259 n = 2;
260 PetscCall(PetscOptionsRealArray("-du", "Perturbation in initial solution", "", du, &n, &flg2));
261 u[0] += du[0];
262 u[1] += du[1];
263 PetscCall(VecRestoreArray(U, &u));
264 if (flg1 || flg2) {
265 ctx.tf = -1;
266 ctx.tcl = -1;
267 }
268 }
269 PetscOptionsEnd();
270
271 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
272 Create timestepping solver context
273 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
274 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
275 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
276 PetscCall(TSSetEquationType(ts, TS_EQ_ODE_EXPLICIT)); /* less Jacobian evaluations when adjoint BEuler is used, otherwise no effect */
277 PetscCall(TSSetType(ts, TSRK));
278 PetscCall(TSSetRHSFunction(ts, NULL, (TSRHSFunctionFn *)RHSFunction, &ctx));
279 PetscCall(TSSetRHSJacobian(ts, A, A, (TSRHSJacobianFn *)RHSJacobian, &ctx));
280 PetscCall(TSCreateQuadratureTS(ts, PETSC_TRUE, &quadts));
281 PetscCall(TSSetRHSFunction(quadts, NULL, (TSRHSFunctionFn *)CostIntegrand, &ctx));
282 PetscCall(TSSetRHSJacobian(quadts, DRDU, DRDU, (TSRHSJacobianFn *)DRDUJacobianTranspose, &ctx));
283 PetscCall(TSSetRHSJacobianP(quadts, DRDP, (TSRHSJacobianPFn *)DRDPJacobianTranspose, &ctx));
284 PetscCall(TSSetCostGradients(ts, 1, lambda, mu));
285 PetscCall(TSSetRHSJacobianP(ts, Jacp, RHSJacobianP, &ctx));
286
287 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
288 Set initial conditions
289 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
290 PetscCall(TSSetSolution(ts, U));
291
292 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
293 Save trajectory of solution so that TSAdjointSolve() may be used
294 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
295 PetscCall(TSSetSaveTrajectory(ts));
296
297 PetscCall(MatCreateVecs(A, &lambda[0], NULL));
298 /* Set initial conditions for the adjoint integration */
299 PetscCall(VecGetArray(lambda[0], &y_ptr));
300 y_ptr[0] = 0.0;
301 y_ptr[1] = 0.0;
302 PetscCall(VecRestoreArray(lambda[0], &y_ptr));
303
304 PetscCall(MatCreateVecs(Jacp, &mu[0], NULL));
305 PetscCall(VecGetArray(mu[0], &x_ptr));
306 x_ptr[0] = -1.0;
307 PetscCall(VecRestoreArray(mu[0], &x_ptr));
308
309 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
310 Set solver options
311 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
312 PetscCall(TSSetMaxTime(ts, 10.0));
313 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
314 PetscCall(TSSetTimeStep(ts, .01));
315 PetscCall(TSSetFromOptions(ts));
316
317 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
318 Solve nonlinear system
319 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
320 if (ensemble) {
321 for (du[1] = -2.5; du[1] <= .01; du[1] += .1) {
322 PetscCall(VecGetArray(U, &u));
323 u[0] = PetscAsinScalar(ctx.Pm / ctx.Pmax);
324 u[1] = ctx.omega_s;
325 u[0] += du[0];
326 u[1] += du[1];
327 PetscCall(VecRestoreArray(U, &u));
328 PetscCall(TSSetTimeStep(ts, .01));
329 PetscCall(TSSolve(ts, U));
330 }
331 } else {
332 PetscCall(TSSolve(ts, U));
333 }
334 PetscCall(VecView(U, PETSC_VIEWER_STDOUT_WORLD));
335 PetscCall(TSGetSolveTime(ts, &ftime));
336 PetscCall(TSGetStepNumber(ts, &steps));
337
338 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
339 Adjoint model starts here
340 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
341 /* Set initial conditions for the adjoint integration */
342 PetscCall(VecGetArray(lambda[0], &y_ptr));
343 y_ptr[0] = 0.0;
344 y_ptr[1] = 0.0;
345 PetscCall(VecRestoreArray(lambda[0], &y_ptr));
346
347 PetscCall(VecGetArray(mu[0], &x_ptr));
348 x_ptr[0] = -1.0;
349 PetscCall(VecRestoreArray(mu[0], &x_ptr));
350
351 PetscCall(TSAdjointSolve(ts));
352
353 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt initial conditions: d[Psi(tf)]/d[phi0] d[Psi(tf)]/d[omega0]\n"));
354 PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
355 PetscCall(VecView(mu[0], PETSC_VIEWER_STDOUT_WORLD));
356 PetscCall(TSGetCostIntegral(ts, &q));
357 PetscCall(VecView(q, PETSC_VIEWER_STDOUT_WORLD));
358 PetscCall(VecGetArray(q, &x_ptr));
359 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n cost function=%g\n", (double)(x_ptr[0] - ctx.Pm)));
360 PetscCall(VecRestoreArray(q, &x_ptr));
361
362 PetscCall(ComputeSensiP(lambda[0], mu[0], &ctx));
363
364 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
365 Free work space. All PETSc objects should be destroyed when they are no longer needed.
366 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
367 PetscCall(MatDestroy(&A));
368 PetscCall(MatDestroy(&Jacp));
369 PetscCall(MatDestroy(&DRDU));
370 PetscCall(MatDestroy(&DRDP));
371 PetscCall(VecDestroy(&U));
372 PetscCall(VecDestroy(&lambda[0]));
373 PetscCall(VecDestroy(&mu[0]));
374 PetscCall(TSDestroy(&ts));
375 PetscCall(PetscFinalize());
376 return 0;
377 }
378
379 /*TEST
380
381 build:
382 requires: !complex
383
384 test:
385 args: -viewer_binary_skip_info -ts_adapt_type none
386
387 TEST*/
388