1 static char help[] = "Sensitivity analysis applied in power grid stability analysis of WECC 9 bus system.\n\
2 This example is based on the 9-bus (node) example given in the book Power\n\
3 Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\
4 The power grid in this example consists of 9 buses (nodes), 3 generators,\n\
5 3 loads, and 9 transmission lines. The network equations are written\n\
6 in current balance form using rectangular coordinates.\n\n";
7
8 /*
9 The equations for the stability analysis are described by the DAE. See ex9bus.c for details.
10 The system has 'jumps' due to faults, thus the time interval is split into multiple sections, and TSSolve is called for each of them. But TSAdjointSolve only needs to be called once since the whole trajectory has been saved in the forward run.
11 The code computes the sensitivity of a final state w.r.t. initial conditions.
12 */
13
14 #include <petscts.h>
15 #include <petscdm.h>
16 #include <petscdmda.h>
17 #include <petscdmcomposite.h>
18
19 #define freq 60
20 #define w_s (2 * PETSC_PI * freq)
21
22 /* Sizes and indices */
23 const PetscInt nbus = 9; /* Number of network buses */
24 const PetscInt ngen = 3; /* Number of generators */
25 const PetscInt nload = 3; /* Number of loads */
26 const PetscInt gbus[3] = {0, 1, 2}; /* Buses at which generators are incident */
27 const PetscInt lbus[3] = {4, 5, 7}; /* Buses at which loads are incident */
28
29 /* Generator real and reactive powers (found via loadflow) */
30 const PetscScalar PG[3] = {0.716786142395021, 1.630000000000000, 0.850000000000000};
31 const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
32 /* Generator constants */
33 const PetscScalar H[3] = {23.64, 6.4, 3.01}; /* Inertia constant */
34 const PetscScalar Rs[3] = {0.0, 0.0, 0.0}; /* Stator Resistance */
35 const PetscScalar Xd[3] = {0.146, 0.8958, 1.3125}; /* d-axis reactance */
36 const PetscScalar Xdp[3] = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
37 const PetscScalar Xq[3] = {0.4360, 0.8645, 1.2578}; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
38 const PetscScalar Xqp[3] = {0.0969, 0.1969, 0.25}; /* q-axis transient reactance */
39 const PetscScalar Td0p[3] = {8.96, 6.0, 5.89}; /* d-axis open circuit time constant */
40 const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6}; /* q-axis open circuit time constant */
41 PetscScalar M[3]; /* M = 2*H/w_s */
42 PetscScalar D[3]; /* D = 0.1*M */
43
44 PetscScalar TM[3]; /* Mechanical Torque */
45 /* Exciter system constants */
46 const PetscScalar KA[3] = {20.0, 20.0, 20.0}; /* Voltage regulartor gain constant */
47 const PetscScalar TA[3] = {0.2, 0.2, 0.2}; /* Voltage regulator time constant */
48 const PetscScalar KE[3] = {1.0, 1.0, 1.0}; /* Exciter gain constant */
49 const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
50 const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
51 const PetscScalar TF[3] = {0.35, 0.35, 0.35}; /* Feedback stabilizer time constant */
52 const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
53 const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */
54
55 PetscScalar Vref[3];
56 /* Load constants
57 We use a composite load model that describes the load and reactive powers at each time instant as follows
58 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
59 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
60 where
61 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
62 ld_alphap,ld_alphap - Percentage contribution (weights) or loads
63 P_D0 - Real power load
64 Q_D0 - Reactive power load
65 V_m(t) - Voltage magnitude at time t
66 V_m0 - Voltage magnitude at t = 0
67 ld_betap, ld_betaq - exponents describing the load model for real and reactive part
68
69 Note: All loads have the same characteristic currently.
70 */
71 const PetscScalar PD0[3] = {1.25, 0.9, 1.0};
72 const PetscScalar QD0[3] = {0.5, 0.3, 0.35};
73 const PetscInt ld_nsegsp[3] = {3, 3, 3};
74 const PetscScalar ld_alphap[3] = {1.0, 0.0, 0.0};
75 const PetscScalar ld_betap[3] = {2.0, 1.0, 0.0};
76 const PetscInt ld_nsegsq[3] = {3, 3, 3};
77 const PetscScalar ld_alphaq[3] = {1.0, 0.0, 0.0};
78 const PetscScalar ld_betaq[3] = {2.0, 1.0, 0.0};
79
80 typedef struct {
81 DM dmgen, dmnet; /* DMs to manage generator and network subsystem */
82 DM dmpgrid; /* Composite DM to manage the entire power grid */
83 Mat Ybus; /* Network admittance matrix */
84 Vec V0; /* Initial voltage vector (Power flow solution) */
85 PetscReal tfaulton, tfaultoff; /* Fault on and off times */
86 PetscInt faultbus; /* Fault bus */
87 PetscScalar Rfault;
88 PetscReal t0, tmax;
89 PetscInt neqs_gen, neqs_net, neqs_pgrid;
90 PetscBool alg_flg;
91 PetscReal t;
92 IS is_diff; /* indices for differential equations */
93 IS is_alg; /* indices for algebraic equations */
94 } Userctx;
95
96 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar * Fr,PetscScalar * Fi)97 PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
98 {
99 PetscFunctionBegin;
100 *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
101 *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
102 PetscFunctionReturn(PETSC_SUCCESS);
103 }
104
105 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar * Fd,PetscScalar * Fq)106 PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
107 {
108 PetscFunctionBegin;
109 *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
110 *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
111 PetscFunctionReturn(PETSC_SUCCESS);
112 }
113
SetInitialGuess(Vec X,Userctx * user)114 PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
115 {
116 Vec Xgen, Xnet;
117 PetscScalar *xgen, *xnet;
118 PetscInt i, idx = 0;
119 PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
120 PetscScalar Eqp, Edp, delta;
121 PetscScalar Efd, RF, VR; /* Exciter variables */
122 PetscScalar Id, Iq; /* Generator dq axis currents */
123 PetscScalar theta, Vd, Vq, SE;
124
125 PetscFunctionBegin;
126 M[0] = 2 * H[0] / w_s;
127 M[1] = 2 * H[1] / w_s;
128 M[2] = 2 * H[2] / w_s;
129 D[0] = 0.1 * M[0];
130 D[1] = 0.1 * M[1];
131 D[2] = 0.1 * M[2];
132
133 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
134
135 /* Network subsystem initialization */
136 PetscCall(VecCopy(user->V0, Xnet));
137
138 /* Generator subsystem initialization */
139 PetscCall(VecGetArray(Xgen, &xgen));
140 PetscCall(VecGetArray(Xnet, &xnet));
141
142 for (i = 0; i < ngen; i++) {
143 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
144 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
145 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
146 Vm2 = Vm * Vm;
147 IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
148 IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;
149
150 delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
151
152 theta = PETSC_PI / 2.0 - delta;
153
154 Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
155 Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
156
157 Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
158 Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
159
160 Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
161 Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
162
163 TM[i] = PG[i];
164
165 /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
166 xgen[idx] = Eqp;
167 xgen[idx + 1] = Edp;
168 xgen[idx + 2] = delta;
169 xgen[idx + 3] = w_s;
170
171 idx = idx + 4;
172
173 xgen[idx] = Id;
174 xgen[idx + 1] = Iq;
175
176 idx = idx + 2;
177
178 /* Exciter */
179 Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
180 SE = k1[i] * PetscExpScalar(k2[i] * Efd);
181 VR = KE[i] * Efd + SE;
182 RF = KF[i] * Efd / TF[i];
183
184 xgen[idx] = Efd;
185 xgen[idx + 1] = RF;
186 xgen[idx + 2] = VR;
187
188 Vref[i] = Vm + (VR / KA[i]);
189
190 idx = idx + 3;
191 }
192
193 PetscCall(VecRestoreArray(Xgen, &xgen));
194 PetscCall(VecRestoreArray(Xnet, &xnet));
195
196 /* PetscCall(VecView(Xgen,0)); */
197 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet));
198 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
199 PetscFunctionReturn(PETSC_SUCCESS);
200 }
201
202 /* Computes F = [f(x,y);g(x,y)] */
ResidualFunction(SNES snes,Vec X,Vec F,Userctx * user)203 PetscErrorCode ResidualFunction(SNES snes, Vec X, Vec F, Userctx *user)
204 {
205 Vec Xgen, Xnet, Fgen, Fnet;
206 PetscScalar *xgen, *xnet, *fgen, *fnet;
207 PetscInt i, idx = 0;
208 PetscScalar Vr, Vi, Vm, Vm2;
209 PetscScalar Eqp, Edp, delta, w; /* Generator variables */
210 PetscScalar Efd, RF, VR; /* Exciter variables */
211 PetscScalar Id, Iq; /* Generator dq axis currents */
212 PetscScalar Vd, Vq, SE;
213 PetscScalar IGr, IGi, IDr, IDi;
214 PetscScalar Zdq_inv[4], det;
215 PetscScalar PD, QD, Vm0, *v0;
216 PetscInt k;
217
218 PetscFunctionBegin;
219 PetscCall(VecZeroEntries(F));
220 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
221 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet));
222 PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
223 PetscCall(DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet));
224
225 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
226 The generator current injection, IG, and load current injection, ID are added later
227 */
228 /* Note that the values in Ybus are stored assuming the imaginary current balance
229 equation is ordered first followed by real current balance equation for each bus.
230 Thus imaginary current contribution goes in location 2*i, and
231 real current contribution in 2*i+1
232 */
233 PetscCall(MatMult(user->Ybus, Xnet, Fnet));
234
235 PetscCall(VecGetArray(Xgen, &xgen));
236 PetscCall(VecGetArray(Xnet, &xnet));
237 PetscCall(VecGetArray(Fgen, &fgen));
238 PetscCall(VecGetArray(Fnet, &fnet));
239
240 /* Generator subsystem */
241 for (i = 0; i < ngen; i++) {
242 Eqp = xgen[idx];
243 Edp = xgen[idx + 1];
244 delta = xgen[idx + 2];
245 w = xgen[idx + 3];
246 Id = xgen[idx + 4];
247 Iq = xgen[idx + 5];
248 Efd = xgen[idx + 6];
249 RF = xgen[idx + 7];
250 VR = xgen[idx + 8];
251
252 /* Generator differential equations */
253 fgen[idx] = (Eqp + (Xd[i] - Xdp[i]) * Id - Efd) / Td0p[i];
254 fgen[idx + 1] = (Edp - (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
255 fgen[idx + 2] = -w + w_s;
256 fgen[idx + 3] = (-TM[i] + Edp * Id + Eqp * Iq + (Xqp[i] - Xdp[i]) * Id * Iq + D[i] * (w - w_s)) / M[i];
257
258 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
259 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
260
261 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
262 /* Algebraic equations for stator currents */
263
264 det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
265
266 Zdq_inv[0] = Rs[i] / det;
267 Zdq_inv[1] = Xqp[i] / det;
268 Zdq_inv[2] = -Xdp[i] / det;
269 Zdq_inv[3] = Rs[i] / det;
270
271 fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
272 fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
273
274 /* Add generator current injection to network */
275 PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
276
277 fnet[2 * gbus[i]] -= IGi;
278 fnet[2 * gbus[i] + 1] -= IGr;
279
280 Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
281
282 SE = k1[i] * PetscExpScalar(k2[i] * Efd);
283
284 /* Exciter differential equations */
285 fgen[idx + 6] = (KE[i] * Efd + SE - VR) / TE[i];
286 fgen[idx + 7] = (RF - KF[i] * Efd / TF[i]) / TF[i];
287 fgen[idx + 8] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
288
289 idx = idx + 9;
290 }
291
292 PetscCall(VecGetArray(user->V0, &v0));
293 for (i = 0; i < nload; i++) {
294 Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
295 Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
296 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
297 Vm2 = Vm * Vm;
298 Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
299 PD = QD = 0.0;
300 for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
301 for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
302
303 /* Load currents */
304 IDr = (PD * Vr + QD * Vi) / Vm2;
305 IDi = (-QD * Vr + PD * Vi) / Vm2;
306
307 fnet[2 * lbus[i]] += IDi;
308 fnet[2 * lbus[i] + 1] += IDr;
309 }
310 PetscCall(VecRestoreArray(user->V0, &v0));
311
312 PetscCall(VecRestoreArray(Xgen, &xgen));
313 PetscCall(VecRestoreArray(Xnet, &xnet));
314 PetscCall(VecRestoreArray(Fgen, &fgen));
315 PetscCall(VecRestoreArray(Fnet, &fnet));
316
317 PetscCall(DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet));
318 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
319 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet));
320 PetscFunctionReturn(PETSC_SUCCESS);
321 }
322
323 /* \dot{x} - f(x,y)
324 g(x,y) = 0
325 */
IFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx * user)326 PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
327 {
328 SNES snes;
329 PetscScalar *f;
330 const PetscScalar *xdot;
331 PetscInt i;
332
333 PetscFunctionBegin;
334 user->t = t;
335
336 PetscCall(TSGetSNES(ts, &snes));
337 PetscCall(ResidualFunction(snes, X, F, user));
338 PetscCall(VecGetArray(F, &f));
339 PetscCall(VecGetArrayRead(Xdot, &xdot));
340 for (i = 0; i < ngen; i++) {
341 f[9 * i] += xdot[9 * i];
342 f[9 * i + 1] += xdot[9 * i + 1];
343 f[9 * i + 2] += xdot[9 * i + 2];
344 f[9 * i + 3] += xdot[9 * i + 3];
345 f[9 * i + 6] += xdot[9 * i + 6];
346 f[9 * i + 7] += xdot[9 * i + 7];
347 f[9 * i + 8] += xdot[9 * i + 8];
348 }
349 PetscCall(VecRestoreArray(F, &f));
350 PetscCall(VecRestoreArrayRead(Xdot, &xdot));
351 PetscFunctionReturn(PETSC_SUCCESS);
352 }
353
354 /* This function is used for solving the algebraic system only during fault on and
355 off times. It computes the entire F and then zeros out the part corresponding to
356 differential equations
357 F = [0;g(y)];
358 */
AlgFunction(SNES snes,Vec X,Vec F,PetscCtx ctx)359 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, PetscCtx ctx)
360 {
361 Userctx *user = (Userctx *)ctx;
362 PetscScalar *f;
363 PetscInt i;
364
365 PetscFunctionBegin;
366 PetscCall(ResidualFunction(snes, X, F, user));
367 PetscCall(VecGetArray(F, &f));
368 for (i = 0; i < ngen; i++) {
369 f[9 * i] = 0;
370 f[9 * i + 1] = 0;
371 f[9 * i + 2] = 0;
372 f[9 * i + 3] = 0;
373 f[9 * i + 6] = 0;
374 f[9 * i + 7] = 0;
375 f[9 * i + 8] = 0;
376 }
377 PetscCall(VecRestoreArray(F, &f));
378 PetscFunctionReturn(PETSC_SUCCESS);
379 }
380
PreallocateJacobian(Mat J,Userctx * user)381 PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
382 {
383 PetscInt *d_nnz;
384 PetscInt i, idx = 0, start = 0;
385 PetscInt ncols;
386
387 PetscFunctionBegin;
388 PetscCall(PetscMalloc1(user->neqs_pgrid, &d_nnz));
389 for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
390 /* Generator subsystem */
391 for (i = 0; i < ngen; i++) {
392 d_nnz[idx] += 3;
393 d_nnz[idx + 1] += 2;
394 d_nnz[idx + 2] += 2;
395 d_nnz[idx + 3] += 5;
396 d_nnz[idx + 4] += 6;
397 d_nnz[idx + 5] += 6;
398
399 d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
400 d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
401
402 d_nnz[idx + 6] += 2;
403 d_nnz[idx + 7] += 2;
404 d_nnz[idx + 8] += 5;
405
406 idx = idx + 9;
407 }
408
409 start = user->neqs_gen;
410
411 for (i = 0; i < nbus; i++) {
412 PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
413 d_nnz[start + 2 * i] += ncols;
414 d_nnz[start + 2 * i + 1] += ncols;
415 PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL));
416 }
417
418 PetscCall(MatSeqAIJSetPreallocation(J, 0, d_nnz));
419
420 PetscCall(PetscFree(d_nnz));
421 PetscFunctionReturn(PETSC_SUCCESS);
422 }
423
424 /*
425 J = [-df_dx, -df_dy
426 dg_dx, dg_dy]
427 */
ResidualJacobian(SNES snes,Vec X,Mat J,Mat B,PetscCtx ctx)428 PetscErrorCode ResidualJacobian(SNES snes, Vec X, Mat J, Mat B, PetscCtx ctx)
429 {
430 Userctx *user = (Userctx *)ctx;
431 Vec Xgen, Xnet;
432 PetscScalar *xgen, *xnet;
433 PetscInt i, idx = 0;
434 PetscScalar Vr, Vi, Vm, Vm2;
435 PetscScalar Eqp, Edp, delta; /* Generator variables */
436 PetscScalar Efd; /* Exciter variables */
437 PetscScalar Id, Iq; /* Generator dq axis currents */
438 PetscScalar Vd, Vq;
439 PetscScalar val[10];
440 PetscInt row[2], col[10];
441 PetscInt net_start = user->neqs_gen;
442 PetscScalar Zdq_inv[4], det;
443 PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
444 PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
445 PetscScalar dSE_dEfd;
446 PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
447 PetscInt ncols;
448 const PetscInt *cols;
449 const PetscScalar *yvals;
450 PetscInt k;
451 PetscScalar PD, QD, Vm0, *v0, Vm4;
452 PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
453 PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
454
455 PetscFunctionBegin;
456 PetscCall(MatZeroEntries(B));
457 PetscCall(DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet));
458 PetscCall(DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet));
459
460 PetscCall(VecGetArray(Xgen, &xgen));
461 PetscCall(VecGetArray(Xnet, &xnet));
462
463 /* Generator subsystem */
464 for (i = 0; i < ngen; i++) {
465 Eqp = xgen[idx];
466 Edp = xgen[idx + 1];
467 delta = xgen[idx + 2];
468 Id = xgen[idx + 4];
469 Iq = xgen[idx + 5];
470 Efd = xgen[idx + 6];
471
472 /* fgen[idx] = (Eqp + (Xd[i] - Xdp[i])*Id - Efd)/Td0p[i]; */
473 row[0] = idx;
474 col[0] = idx;
475 col[1] = idx + 4;
476 col[2] = idx + 6;
477 val[0] = 1 / Td0p[i];
478 val[1] = (Xd[i] - Xdp[i]) / Td0p[i];
479 val[2] = -1 / Td0p[i];
480
481 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
482
483 /* fgen[idx+1] = (Edp - (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
484 row[0] = idx + 1;
485 col[0] = idx + 1;
486 col[1] = idx + 5;
487 val[0] = 1 / Tq0p[i];
488 val[1] = -(Xq[i] - Xqp[i]) / Tq0p[i];
489 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
490
491 /* fgen[idx+2] = - w + w_s; */
492 row[0] = idx + 2;
493 col[0] = idx + 2;
494 col[1] = idx + 3;
495 val[0] = 0;
496 val[1] = -1;
497 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
498
499 /* fgen[idx+3] = (-TM[i] + Edp*Id + Eqp*Iq + (Xqp[i] - Xdp[i])*Id*Iq + D[i]*(w - w_s))/M[i]; */
500 row[0] = idx + 3;
501 col[0] = idx;
502 col[1] = idx + 1;
503 col[2] = idx + 3;
504 col[3] = idx + 4;
505 col[4] = idx + 5;
506 val[0] = Iq / M[i];
507 val[1] = Id / M[i];
508 val[2] = D[i] / M[i];
509 val[3] = (Edp + (Xqp[i] - Xdp[i]) * Iq) / M[i];
510 val[4] = (Eqp + (Xqp[i] - Xdp[i]) * Id) / M[i];
511 PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
512
513 Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
514 Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
515 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
516
517 det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
518
519 Zdq_inv[0] = Rs[i] / det;
520 Zdq_inv[1] = Xqp[i] / det;
521 Zdq_inv[2] = -Xdp[i] / det;
522 Zdq_inv[3] = Rs[i] / det;
523
524 dVd_dVr = PetscSinScalar(delta);
525 dVd_dVi = -PetscCosScalar(delta);
526 dVq_dVr = PetscCosScalar(delta);
527 dVq_dVi = PetscSinScalar(delta);
528 dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
529 dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
530
531 /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
532 row[0] = idx + 4;
533 col[0] = idx;
534 col[1] = idx + 1;
535 col[2] = idx + 2;
536 val[0] = -Zdq_inv[1];
537 val[1] = -Zdq_inv[0];
538 val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
539 col[3] = idx + 4;
540 col[4] = net_start + 2 * gbus[i];
541 col[5] = net_start + 2 * gbus[i] + 1;
542 val[3] = 1;
543 val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
544 val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
545 PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
546
547 /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
548 row[0] = idx + 5;
549 col[0] = idx;
550 col[1] = idx + 1;
551 col[2] = idx + 2;
552 val[0] = -Zdq_inv[3];
553 val[1] = -Zdq_inv[2];
554 val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
555 col[3] = idx + 5;
556 col[4] = net_start + 2 * gbus[i];
557 col[5] = net_start + 2 * gbus[i] + 1;
558 val[3] = 1;
559 val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
560 val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
561 PetscCall(MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES));
562
563 dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
564 dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
565 dIGr_dId = PetscSinScalar(delta);
566 dIGr_dIq = PetscCosScalar(delta);
567 dIGi_dId = -PetscCosScalar(delta);
568 dIGi_dIq = PetscSinScalar(delta);
569
570 /* fnet[2*gbus[i]] -= IGi; */
571 row[0] = net_start + 2 * gbus[i];
572 col[0] = idx + 2;
573 col[1] = idx + 4;
574 col[2] = idx + 5;
575 val[0] = -dIGi_ddelta;
576 val[1] = -dIGi_dId;
577 val[2] = -dIGi_dIq;
578 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
579
580 /* fnet[2*gbus[i]+1] -= IGr; */
581 row[0] = net_start + 2 * gbus[i] + 1;
582 col[0] = idx + 2;
583 col[1] = idx + 4;
584 col[2] = idx + 5;
585 val[0] = -dIGr_ddelta;
586 val[1] = -dIGr_dId;
587 val[2] = -dIGr_dIq;
588 PetscCall(MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES));
589
590 Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
591
592 /* fgen[idx+6] = (KE[i]*Efd + SE - VR)/TE[i]; */
593 /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
594 dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
595
596 row[0] = idx + 6;
597 col[0] = idx + 6;
598 col[1] = idx + 8;
599 val[0] = (KE[i] + dSE_dEfd) / TE[i];
600 val[1] = -1 / TE[i];
601 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
602
603 /* Exciter differential equations */
604
605 /* fgen[idx+7] = (RF - KF[i]*Efd/TF[i])/TF[i]; */
606 row[0] = idx + 7;
607 col[0] = idx + 6;
608 col[1] = idx + 7;
609 val[0] = (-KF[i] / TF[i]) / TF[i];
610 val[1] = 1 / TF[i];
611 PetscCall(MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES));
612
613 /* fgen[idx+8] = (VR - KA[i]*RF + KA[i]*KF[i]*Efd/TF[i] - KA[i]*(Vref[i] - Vm))/TA[i]; */
614 /* Vm = (Vd^2 + Vq^2)^0.5; */
615 dVm_dVd = Vd / Vm;
616 dVm_dVq = Vq / Vm;
617 dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
618 dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
619 row[0] = idx + 8;
620 col[0] = idx + 6;
621 col[1] = idx + 7;
622 col[2] = idx + 8;
623 val[0] = (KA[i] * KF[i] / TF[i]) / TA[i];
624 val[1] = -KA[i] / TA[i];
625 val[2] = 1 / TA[i];
626 col[3] = net_start + 2 * gbus[i];
627 col[4] = net_start + 2 * gbus[i] + 1;
628 val[3] = KA[i] * dVm_dVr / TA[i];
629 val[4] = KA[i] * dVm_dVi / TA[i];
630 PetscCall(MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES));
631 idx = idx + 9;
632 }
633
634 for (i = 0; i < nbus; i++) {
635 PetscCall(MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
636 row[0] = net_start + 2 * i;
637 for (k = 0; k < ncols; k++) {
638 col[k] = net_start + cols[k];
639 val[k] = yvals[k];
640 }
641 PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
642 PetscCall(MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals));
643
644 PetscCall(MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
645 row[0] = net_start + 2 * i + 1;
646 for (k = 0; k < ncols; k++) {
647 col[k] = net_start + cols[k];
648 val[k] = yvals[k];
649 }
650 PetscCall(MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES));
651 PetscCall(MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals));
652 }
653
654 PetscCall(MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY));
655 PetscCall(MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY));
656
657 PetscCall(VecGetArray(user->V0, &v0));
658 for (i = 0; i < nload; i++) {
659 Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
660 Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
661 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
662 Vm2 = Vm * Vm;
663 Vm4 = Vm2 * Vm2;
664 Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
665 PD = QD = 0.0;
666 dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
667 for (k = 0; k < ld_nsegsp[i]; k++) {
668 PD += ld_alphap[k] * PD0[i] * PetscPowScalar(Vm / Vm0, ld_betap[k]);
669 dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vr * PetscPowScalar(Vm, ld_betap[k] - 2);
670 dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar(1 / Vm0, ld_betap[k]) * Vi * PetscPowScalar(Vm, ld_betap[k] - 2);
671 }
672 for (k = 0; k < ld_nsegsq[i]; k++) {
673 QD += ld_alphaq[k] * QD0[i] * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
674 dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vr * PetscPowScalar(Vm, ld_betaq[k] - 2);
675 dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar(1 / Vm0, ld_betaq[k]) * Vi * PetscPowScalar(Vm, ld_betaq[k] - 2);
676 }
677
678 /* IDr = (PD*Vr + QD*Vi)/Vm2; */
679 /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
680
681 dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
682 dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
683
684 dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
685 dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
686
687 /* fnet[2*lbus[i]] += IDi; */
688 row[0] = net_start + 2 * lbus[i];
689 col[0] = net_start + 2 * lbus[i];
690 col[1] = net_start + 2 * lbus[i] + 1;
691 val[0] = dIDi_dVr;
692 val[1] = dIDi_dVi;
693 PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
694 /* fnet[2*lbus[i]+1] += IDr; */
695 row[0] = net_start + 2 * lbus[i] + 1;
696 col[0] = net_start + 2 * lbus[i];
697 col[1] = net_start + 2 * lbus[i] + 1;
698 val[0] = dIDr_dVr;
699 val[1] = dIDr_dVi;
700 PetscCall(MatSetValues(J, 1, row, 2, col, val, ADD_VALUES));
701 }
702 PetscCall(VecRestoreArray(user->V0, &v0));
703
704 PetscCall(VecRestoreArray(Xgen, &xgen));
705 PetscCall(VecRestoreArray(Xnet, &xnet));
706
707 PetscCall(DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet));
708
709 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
710 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
711 PetscFunctionReturn(PETSC_SUCCESS);
712 }
713
714 /*
715 J = [I, 0
716 dg_dx, dg_dy]
717 */
AlgJacobian(SNES snes,Vec X,Mat A,Mat B,PetscCtx ctx)718 PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, PetscCtx ctx)
719 {
720 Userctx *user = (Userctx *)ctx;
721
722 PetscFunctionBegin;
723 PetscCall(ResidualJacobian(snes, X, A, B, ctx));
724 PetscCall(MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE));
725 PetscCall(MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL));
726 PetscFunctionReturn(PETSC_SUCCESS);
727 }
728
729 /*
730 J = [a*I-df_dx, -df_dy
731 dg_dx, dg_dy]
732 */
733
IJacobian(TS ts,PetscReal t,Vec X,Vec Xdot,PetscReal a,Mat A,Mat B,Userctx * user)734 PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
735 {
736 SNES snes;
737 PetscScalar atmp = (PetscScalar)a;
738 PetscInt i, row;
739
740 PetscFunctionBegin;
741 user->t = t;
742
743 PetscCall(TSGetSNES(ts, &snes));
744 PetscCall(ResidualJacobian(snes, X, A, B, user));
745 for (i = 0; i < ngen; i++) {
746 row = 9 * i;
747 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
748 row = 9 * i + 1;
749 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
750 row = 9 * i + 2;
751 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
752 row = 9 * i + 3;
753 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
754 row = 9 * i + 6;
755 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
756 row = 9 * i + 7;
757 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
758 row = 9 * i + 8;
759 PetscCall(MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES));
760 }
761 PetscCall(MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY));
762 PetscCall(MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY));
763 PetscFunctionReturn(PETSC_SUCCESS);
764 }
765
main(int argc,char ** argv)766 int main(int argc, char **argv)
767 {
768 TS ts;
769 SNES snes_alg;
770 PetscMPIInt size;
771 Userctx user;
772 PetscViewer Xview, Ybusview;
773 Vec X;
774 Mat J;
775 PetscInt i;
776 /* sensitivity context */
777 PetscScalar *y_ptr;
778 Vec lambda[1];
779 PetscInt *idx2;
780 Vec Xdot;
781 Vec F_alg;
782 PetscInt row_loc, col_loc;
783 PetscScalar val;
784
785 PetscFunctionBeginUser;
786 PetscCall(PetscInitialize(&argc, &argv, "petscoptions", help));
787 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
788 PetscCheck(size == 1, PETSC_COMM_WORLD, PETSC_ERR_WRONG_MPI_SIZE, "Only for sequential runs");
789
790 user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
791 user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
792 user.neqs_pgrid = user.neqs_gen + user.neqs_net;
793
794 /* Create indices for differential and algebraic equations */
795 PetscCall(PetscMalloc1(7 * ngen, &idx2));
796 for (i = 0; i < ngen; i++) {
797 idx2[7 * i] = 9 * i;
798 idx2[7 * i + 1] = 9 * i + 1;
799 idx2[7 * i + 2] = 9 * i + 2;
800 idx2[7 * i + 3] = 9 * i + 3;
801 idx2[7 * i + 4] = 9 * i + 6;
802 idx2[7 * i + 5] = 9 * i + 7;
803 idx2[7 * i + 6] = 9 * i + 8;
804 }
805 PetscCall(ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff));
806 PetscCall(ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg));
807 PetscCall(PetscFree(idx2));
808
809 /* Read initial voltage vector and Ybus */
810 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview));
811 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview));
812
813 PetscCall(VecCreate(PETSC_COMM_WORLD, &user.V0));
814 PetscCall(VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net));
815 PetscCall(VecLoad(user.V0, Xview));
816
817 PetscCall(MatCreate(PETSC_COMM_WORLD, &user.Ybus));
818 PetscCall(MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net));
819 PetscCall(MatSetType(user.Ybus, MATBAIJ));
820 /* PetscCall(MatSetBlockSize(user.Ybus,2)); */
821 PetscCall(MatLoad(user.Ybus, Ybusview));
822
823 /* Set run time options */
824 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
825 {
826 user.tfaulton = 1.0;
827 user.tfaultoff = 1.2;
828 user.Rfault = 0.0001;
829 user.faultbus = 8;
830 PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
831 PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
832 PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
833 user.t0 = 0.0;
834 user.tmax = 5.0;
835 PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
836 PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
837 }
838 PetscOptionsEnd();
839
840 PetscCall(PetscViewerDestroy(&Xview));
841 PetscCall(PetscViewerDestroy(&Ybusview));
842
843 /* Create DMs for generator and network subsystems */
844 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen));
845 PetscCall(DMSetOptionsPrefix(user.dmgen, "dmgen_"));
846 PetscCall(DMSetFromOptions(user.dmgen));
847 PetscCall(DMSetUp(user.dmgen));
848 PetscCall(DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet));
849 PetscCall(DMSetOptionsPrefix(user.dmnet, "dmnet_"));
850 PetscCall(DMSetFromOptions(user.dmnet));
851 PetscCall(DMSetUp(user.dmnet));
852 /* Create a composite DM packer and add the two DMs */
853 PetscCall(DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid));
854 PetscCall(DMSetOptionsPrefix(user.dmpgrid, "pgrid_"));
855 PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmgen));
856 PetscCall(DMCompositeAddDM(user.dmpgrid, user.dmnet));
857
858 PetscCall(DMCreateGlobalVector(user.dmpgrid, &X));
859
860 PetscCall(MatCreate(PETSC_COMM_WORLD, &J));
861 PetscCall(MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid));
862 PetscCall(MatSetFromOptions(J));
863 PetscCall(PreallocateJacobian(J, &user));
864
865 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
866 Create timestepping solver context
867 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
868 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
869 PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
870 PetscCall(TSSetType(ts, TSCN));
871 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)IFunction, &user));
872 PetscCall(TSSetIJacobian(ts, J, J, (TSIJacobianFn *)IJacobian, &user));
873 PetscCall(TSSetApplicationContext(ts, &user));
874
875 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
876 Set initial conditions
877 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
878 PetscCall(SetInitialGuess(X, &user));
879 /* Just to set up the Jacobian structure */
880 PetscCall(VecDuplicate(X, &Xdot));
881 PetscCall(IJacobian(ts, 0.0, X, Xdot, 0.0, J, J, &user));
882 PetscCall(VecDestroy(&Xdot));
883
884 /*
885 Save trajectory of solution so that TSAdjointSolve() may be used
886 */
887 PetscCall(TSSetSaveTrajectory(ts));
888
889 PetscCall(TSSetMaxTime(ts, user.tmax));
890 PetscCall(TSSetTimeStep(ts, 0.01));
891 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));
892 PetscCall(TSSetFromOptions(ts));
893
894 user.alg_flg = PETSC_FALSE;
895 /* Prefault period */
896 PetscCall(TSSolve(ts, X));
897
898 /* Create the nonlinear solver for solving the algebraic system */
899 /* Note that although the algebraic system needs to be solved only for
900 Idq and V, we reuse the entire system including xgen. The xgen
901 variables are held constant by setting their residuals to 0 and
902 putting a 1 on the Jacobian diagonal for xgen rows
903 */
904 PetscCall(VecDuplicate(X, &F_alg));
905 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
906 PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));
907 PetscCall(MatZeroEntries(J));
908 PetscCall(SNESSetJacobian(snes_alg, J, J, AlgJacobian, &user));
909 PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
910 PetscCall(SNESSetFromOptions(snes_alg));
911
912 /* Apply disturbance - resistive fault at user.faultbus */
913 /* This is done by adding shunt conductance to the diagonal location
914 in the Ybus matrix */
915 row_loc = 2 * user.faultbus;
916 col_loc = 2 * user.faultbus + 1; /* Location for G */
917 val = 1 / user.Rfault;
918 PetscCall(MatSetValues(user.Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
919 row_loc = 2 * user.faultbus + 1;
920 col_loc = 2 * user.faultbus; /* Location for G */
921 val = 1 / user.Rfault;
922 PetscCall(MatSetValues(user.Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
923
924 PetscCall(MatAssemblyBegin(user.Ybus, MAT_FINAL_ASSEMBLY));
925 PetscCall(MatAssemblyEnd(user.Ybus, MAT_FINAL_ASSEMBLY));
926
927 user.alg_flg = PETSC_TRUE;
928 /* Solve the algebraic equations */
929 PetscCall(SNESSolve(snes_alg, NULL, X));
930
931 /* Disturbance period */
932 user.alg_flg = PETSC_FALSE;
933 PetscCall(TSSetTime(ts, user.tfaulton));
934 PetscCall(TSSetMaxTime(ts, user.tfaultoff));
935 PetscCall(TSSolve(ts, X));
936
937 /* Remove the fault */
938 row_loc = 2 * user.faultbus;
939 col_loc = 2 * user.faultbus + 1;
940 val = -1 / user.Rfault;
941 PetscCall(MatSetValues(user.Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
942 row_loc = 2 * user.faultbus + 1;
943 col_loc = 2 * user.faultbus;
944 val = -1 / user.Rfault;
945 PetscCall(MatSetValues(user.Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES));
946
947 PetscCall(MatAssemblyBegin(user.Ybus, MAT_FINAL_ASSEMBLY));
948 PetscCall(MatAssemblyEnd(user.Ybus, MAT_FINAL_ASSEMBLY));
949
950 PetscCall(MatZeroEntries(J));
951
952 user.alg_flg = PETSC_TRUE;
953
954 /* Solve the algebraic equations */
955 PetscCall(SNESSolve(snes_alg, NULL, X));
956
957 /* Post-disturbance period */
958 user.alg_flg = PETSC_TRUE;
959 PetscCall(TSSetTime(ts, user.tfaultoff));
960 PetscCall(TSSetMaxTime(ts, user.tmax));
961 PetscCall(TSSolve(ts, X));
962
963 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
964 Adjoint model starts here
965 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
966 PetscCall(TSSetPostStep(ts, NULL));
967 PetscCall(MatCreateVecs(J, &lambda[0], NULL));
968 /* Set initial conditions for the adjoint integration */
969 PetscCall(VecZeroEntries(lambda[0]));
970 PetscCall(VecGetArray(lambda[0], &y_ptr));
971 y_ptr[0] = 1.0;
972 PetscCall(VecRestoreArray(lambda[0], &y_ptr));
973 PetscCall(TSSetCostGradients(ts, 1, lambda, NULL));
974
975 PetscCall(TSAdjointSolve(ts));
976
977 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\n sensitivity wrt initial conditions: \n"));
978 PetscCall(VecView(lambda[0], PETSC_VIEWER_STDOUT_WORLD));
979 PetscCall(VecDestroy(&lambda[0]));
980
981 PetscCall(SNESDestroy(&snes_alg));
982 PetscCall(VecDestroy(&F_alg));
983 PetscCall(MatDestroy(&J));
984 PetscCall(MatDestroy(&user.Ybus));
985 PetscCall(VecDestroy(&X));
986 PetscCall(VecDestroy(&user.V0));
987 PetscCall(DMDestroy(&user.dmgen));
988 PetscCall(DMDestroy(&user.dmnet));
989 PetscCall(DMDestroy(&user.dmpgrid));
990 PetscCall(ISDestroy(&user.is_diff));
991 PetscCall(ISDestroy(&user.is_alg));
992 PetscCall(TSDestroy(&ts));
993 PetscCall(PetscFinalize());
994 return 0;
995 }
996
997 /*TEST
998
999 build:
1000 requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1001
1002 test:
1003 args: -viewer_binary_skip_info
1004 localrunfiles: petscoptions X.bin Ybus.bin
1005
1006 test:
1007 suffix: 2
1008 args: -viewer_binary_skip_info -ts_type beuler
1009 localrunfiles: petscoptions X.bin Ybus.bin
1010
1011 TEST*/
1012