1 static char help[] = "This example uses the same problem set up of ex9busdmnetwork.c. \n\
2 It demonstrates setting and accessing of variables for individual components, instead of \n\
3 the network vertices (as used in ex9busdmnetwork.c). This is especially useful where vertices \n\
4 /edges have multiple-components associated with them and one or more components has physics \n\
5 associated with it. \n\
6 Input parameters include:\n\
7 -nc : number of copies of the base case\n\n";
8
9 /*
10 This example was modified from ex9busdmnetwork.c.
11 */
12
13 #include <petscts.h>
14 #include <petscdmnetwork.h>
15
16 #define FREQ 60
17 #define W_S (2 * PETSC_PI * FREQ)
18 #define NGEN 3 /* No. of generators in the 9 bus system */
19 #define NLOAD 3 /* No. of loads in the 9 bus system */
20 #define NBUS 9 /* No. of buses in the 9 bus system */
21 #define NBRANCH 9 /* No. of branches in the 9 bus system */
22
23 typedef struct {
24 PetscInt id; /* Bus Number or extended bus name*/
25 PetscScalar mbase; /* MVA base of the machine */
26 PetscScalar PG; /* Generator active power output */
27 PetscScalar QG; /* Generator reactive power output */
28
29 /* Generator constants */
30 PetscScalar H; /* Inertia constant */
31 PetscScalar Rs; /* Stator Resistance */
32 PetscScalar Xd; /* d-axis reactance */
33 PetscScalar Xdp; /* d-axis transient reactance */
34 PetscScalar Xq; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
35 PetscScalar Xqp; /* q-axis transient reactance */
36 PetscScalar Td0p; /* d-axis open circuit time constant */
37 PetscScalar Tq0p; /* q-axis open circuit time constant */
38 PetscScalar M; /* M = 2*H/W_S */
39 PetscScalar D; /* D = 0.1*M */
40 PetscScalar TM; /* Mechanical Torque */
41 } Gen;
42
43 typedef struct {
44 /* Exciter system constants */
45 PetscScalar KA; /* Voltage regulator gain constant */
46 PetscScalar TA; /* Voltage regulator time constant */
47 PetscScalar KE; /* Exciter gain constant */
48 PetscScalar TE; /* Exciter time constant */
49 PetscScalar KF; /* Feedback stabilizer gain constant */
50 PetscScalar TF; /* Feedback stabilizer time constant */
51 PetscScalar k1, k2; /* calculating the saturation function SE = k1*exp(k2*Efd) */
52 PetscScalar Vref; /* Voltage regulator voltage setpoint */
53 } Exc;
54
55 typedef struct {
56 PetscInt id; /* node id */
57 PetscInt nofgen; /* Number of generators at the bus*/
58 PetscInt nofload; /* Number of load at the bus*/
59 PetscScalar yff[2]; /* yff[0]= imaginary part of admittance, yff[1]=real part of admittance*/
60 PetscScalar vr; /* Real component of bus voltage */
61 PetscScalar vi; /* Imaginary component of bus voltage */
62 } Bus;
63
64 /* Load constants
65 We use a composite load model that describes the load and reactive powers at each time instant as follows
66 P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
67 Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
68 where
69 id - index of the load
70 ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
71 ld_alphap,ld_alphap - Percentage contribution (weights) or loads
72 P_D0 - Real power load
73 Q_D0 - Reactive power load
74 Vm(t) - Voltage magnitude at time t
75 Vm0 - Voltage magnitude at t = 0
76 ld_betap, ld_betaq - exponents describing the load model for real and reactive part
77
78 Note: All loads have the same characteristic currently.
79 */
80 typedef struct {
81 PetscInt id; /* bus id */
82 PetscInt ld_nsegsp, ld_nsegsq;
83 PetscScalar PD0, QD0;
84 PetscScalar ld_alphap[3]; /* ld_alphap=[1,0,0], an array, not a value, so use *ld_alphap; */
85 PetscScalar ld_betap[3], ld_alphaq[3], ld_betaq[3];
86 } Load;
87
88 typedef struct {
89 PetscInt id; /* node id */
90 PetscScalar yft[2]; /* yft[0]= imaginary part of admittance, yft[1]=real part of admittance*/
91 } Branch;
92
93 typedef struct {
94 PetscReal tfaulton, tfaultoff; /* Fault on and off times */
95 PetscReal t;
96 PetscReal t0, tmax; /* initial time and final time */
97 PetscInt faultbus; /* Fault bus */
98 PetscScalar Rfault; /* Fault resistance (pu) */
99 PetscScalar *ybusfault;
100 PetscBool alg_flg;
101 } Userctx;
102
103 /* Used to read data into the DMNetwork components */
read_data(PetscInt nc,Gen ** pgen,Exc ** pexc,Load ** pload,Bus ** pbus,Branch ** pbranch,PetscInt ** pedgelist)104 PetscErrorCode read_data(PetscInt nc, Gen **pgen, Exc **pexc, Load **pload, Bus **pbus, Branch **pbranch, PetscInt **pedgelist)
105 {
106 PetscInt i, j, row[1], col[2];
107 PetscInt *edgelist;
108 PetscInt nofgen[9] = {1, 1, 1, 0, 0, 0, 0, 0, 0}; /* Buses at which generators are incident */
109 PetscInt nofload[9] = {0, 0, 0, 0, 1, 1, 0, 1, 0}; /* Buses at which loads are incident */
110 const PetscScalar *varr;
111 PetscScalar M[3], D[3];
112 Bus *bus;
113 Branch *branch;
114 Gen *gen;
115 Exc *exc;
116 Load *load;
117 Mat Ybus;
118 Vec V0;
119
120 /*10 parameters*/
121 /* Generator real and reactive powers (found via loadflow) */
122 static const PetscScalar PG[3] = {0.716786142395021, 1.630000000000000, 0.850000000000000};
123 static const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
124
125 /* Generator constants */
126 static const PetscScalar H[3] = {23.64, 6.4, 3.01}; /* Inertia constant */
127 static const PetscScalar Rs[3] = {0.0, 0.0, 0.0}; /* Stator Resistance */
128 static const PetscScalar Xd[3] = {0.146, 0.8958, 1.3125}; /* d-axis reactance */
129 static const PetscScalar Xdp[3] = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
130 static 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 */
131 static const PetscScalar Xqp[3] = {0.0969, 0.1969, 0.25}; /* q-axis transient reactance */
132 static const PetscScalar Td0p[3] = {8.96, 6.0, 5.89}; /* d-axis open circuit time constant */
133 static const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6}; /* q-axis open circuit time constant */
134
135 /* Exciter system constants (8 parameters)*/
136 static const PetscScalar KA[3] = {20.0, 20.0, 20.0}; /* Voltage regulartor gain constant */
137 static const PetscScalar TA[3] = {0.2, 0.2, 0.2}; /* Voltage regulator time constant */
138 static const PetscScalar KE[3] = {1.0, 1.0, 1.0}; /* Exciter gain constant */
139 static const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
140 static const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
141 static const PetscScalar TF[3] = {0.35, 0.35, 0.35}; /* Feedback stabilizer time constant */
142 static const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
143 static const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */
144
145 /* Load constants */
146 static const PetscScalar PD0[3] = {1.25, 0.9, 1.0};
147 static const PetscScalar QD0[3] = {0.5, 0.3, 0.35};
148 static const PetscScalar ld_alphaq[3] = {1, 0, 0};
149 static const PetscScalar ld_betaq[3] = {2, 1, 0};
150 static const PetscScalar ld_betap[3] = {2, 1, 0};
151 static const PetscScalar ld_alphap[3] = {1, 0, 0};
152 PetscInt ld_nsegsp[3] = {3, 3, 3};
153 PetscInt ld_nsegsq[3] = {3, 3, 3};
154 PetscViewer Xview, Ybusview;
155 PetscInt neqs_net, m, n;
156
157 PetscFunctionBeginUser;
158 /* Read V0 and Ybus from files */
159 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "X.bin", FILE_MODE_READ, &Xview));
160 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "Ybus.bin", FILE_MODE_READ, &Ybusview));
161 PetscCall(VecCreate(PETSC_COMM_SELF, &V0));
162 PetscCall(VecLoad(V0, Xview));
163
164 PetscCall(MatCreate(PETSC_COMM_SELF, &Ybus));
165 PetscCall(MatSetType(Ybus, MATBAIJ));
166 PetscCall(MatLoad(Ybus, Ybusview));
167
168 /* Destroy unnecessary stuff */
169 PetscCall(PetscViewerDestroy(&Xview));
170 PetscCall(PetscViewerDestroy(&Ybusview));
171
172 PetscCall(MatGetLocalSize(Ybus, &m, &n));
173 neqs_net = 2 * NBUS; /* # eqs. for network subsystem */
174 PetscCheck(m == neqs_net && n == neqs_net, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix Ybus is in wrong sizes");
175
176 M[0] = 2 * H[0] / W_S;
177 M[1] = 2 * H[1] / W_S;
178 M[2] = 2 * H[2] / W_S;
179 D[0] = 0.1 * M[0];
180 D[1] = 0.1 * M[1];
181 D[2] = 0.1 * M[2];
182
183 /* Allocate memory for bus, generators, exciter, loads and branches */
184 PetscCall(PetscCalloc5(NBUS * nc, &bus, NGEN * nc, &gen, NLOAD * nc, &load, NBRANCH * nc + (nc - 1), &branch, NGEN * nc, &exc));
185
186 PetscCall(VecGetArrayRead(V0, &varr));
187
188 /* read bus data */
189 for (i = 0; i < nc; i++) {
190 for (j = 0; j < NBUS; j++) {
191 bus[i * 9 + j].id = i * 9 + j;
192 bus[i * 9 + j].nofgen = nofgen[j];
193 bus[i * 9 + j].nofload = nofload[j];
194 bus[i * 9 + j].vr = varr[2 * j];
195 bus[i * 9 + j].vi = varr[2 * j + 1];
196 row[0] = 2 * j;
197 col[0] = 2 * j;
198 col[1] = 2 * j + 1;
199 /* real and imaginary part of admittance from Ybus into yff */
200 PetscCall(MatGetValues(Ybus, 1, row, 2, col, bus[i * 9 + j].yff));
201 }
202 }
203
204 /* read generator data */
205 for (i = 0; i < nc; i++) {
206 for (j = 0; j < NGEN; j++) {
207 gen[i * 3 + j].id = i * 3 + j;
208 gen[i * 3 + j].PG = PG[j];
209 gen[i * 3 + j].QG = QG[j];
210 gen[i * 3 + j].H = H[j];
211 gen[i * 3 + j].Rs = Rs[j];
212 gen[i * 3 + j].Xd = Xd[j];
213 gen[i * 3 + j].Xdp = Xdp[j];
214 gen[i * 3 + j].Xq = Xq[j];
215 gen[i * 3 + j].Xqp = Xqp[j];
216 gen[i * 3 + j].Td0p = Td0p[j];
217 gen[i * 3 + j].Tq0p = Tq0p[j];
218 gen[i * 3 + j].M = M[j];
219 gen[i * 3 + j].D = D[j];
220 }
221 }
222
223 for (i = 0; i < nc; i++) {
224 for (j = 0; j < NGEN; j++) {
225 /* exciter system */
226 exc[i * 3 + j].KA = KA[j];
227 exc[i * 3 + j].TA = TA[j];
228 exc[i * 3 + j].KE = KE[j];
229 exc[i * 3 + j].TE = TE[j];
230 exc[i * 3 + j].KF = KF[j];
231 exc[i * 3 + j].TF = TF[j];
232 exc[i * 3 + j].k1 = k1[j];
233 exc[i * 3 + j].k2 = k2[j];
234 }
235 }
236
237 /* read load data */
238 for (i = 0; i < nc; i++) {
239 for (j = 0; j < NLOAD; j++) {
240 load[i * 3 + j].id = i * 3 + j;
241 load[i * 3 + j].PD0 = PD0[j];
242 load[i * 3 + j].QD0 = QD0[j];
243 load[i * 3 + j].ld_nsegsp = ld_nsegsp[j];
244
245 load[i * 3 + j].ld_alphap[0] = ld_alphap[0];
246 load[i * 3 + j].ld_alphap[1] = ld_alphap[1];
247 load[i * 3 + j].ld_alphap[2] = ld_alphap[2];
248
249 load[i * 3 + j].ld_alphaq[0] = ld_alphaq[0];
250 load[i * 3 + j].ld_alphaq[1] = ld_alphaq[1];
251 load[i * 3 + j].ld_alphaq[2] = ld_alphaq[2];
252
253 load[i * 3 + j].ld_betap[0] = ld_betap[0];
254 load[i * 3 + j].ld_betap[1] = ld_betap[1];
255 load[i * 3 + j].ld_betap[2] = ld_betap[2];
256 load[i * 3 + j].ld_nsegsq = ld_nsegsq[j];
257
258 load[i * 3 + j].ld_betaq[0] = ld_betaq[0];
259 load[i * 3 + j].ld_betaq[1] = ld_betaq[1];
260 load[i * 3 + j].ld_betaq[2] = ld_betaq[2];
261 }
262 }
263 PetscCall(PetscCalloc1(2 * NBRANCH * nc + 2 * (nc - 1), &edgelist));
264
265 /* read edgelist */
266 for (i = 0; i < nc; i++) {
267 for (j = 0; j < NBRANCH; j++) {
268 switch (j) {
269 case 0:
270 edgelist[i * 18 + 2 * j] = 0 + 9 * i;
271 edgelist[i * 18 + 2 * j + 1] = 3 + 9 * i;
272 break;
273 case 1:
274 edgelist[i * 18 + 2 * j] = 1 + 9 * i;
275 edgelist[i * 18 + 2 * j + 1] = 6 + 9 * i;
276 break;
277 case 2:
278 edgelist[i * 18 + 2 * j] = 2 + 9 * i;
279 edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
280 break;
281 case 3:
282 edgelist[i * 18 + 2 * j] = 3 + 9 * i;
283 edgelist[i * 18 + 2 * j + 1] = 4 + 9 * i;
284 break;
285 case 4:
286 edgelist[i * 18 + 2 * j] = 3 + 9 * i;
287 edgelist[i * 18 + 2 * j + 1] = 5 + 9 * i;
288 break;
289 case 5:
290 edgelist[i * 18 + 2 * j] = 4 + 9 * i;
291 edgelist[i * 18 + 2 * j + 1] = 6 + 9 * i;
292 break;
293 case 6:
294 edgelist[i * 18 + 2 * j] = 5 + 9 * i;
295 edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
296 break;
297 case 7:
298 edgelist[i * 18 + 2 * j] = 6 + 9 * i;
299 edgelist[i * 18 + 2 * j + 1] = 7 + 9 * i;
300 break;
301 case 8:
302 edgelist[i * 18 + 2 * j] = 7 + 9 * i;
303 edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
304 break;
305 default:
306 break;
307 }
308 }
309 }
310
311 /* for connecting last bus of previous network(9*i-1) to first bus of next network(9*i), the branch admittance=-0.0301407+j17.3611 */
312 for (i = 1; i < nc; i++) {
313 edgelist[18 * nc + 2 * (i - 1)] = 8 + (i - 1) * 9;
314 edgelist[18 * nc + 2 * (i - 1) + 1] = 9 * i;
315
316 /* adding admittances to the off-diagonal elements */
317 branch[9 * nc + (i - 1)].id = 9 * nc + (i - 1);
318 branch[9 * nc + (i - 1)].yft[0] = 17.3611;
319 branch[9 * nc + (i - 1)].yft[1] = -0.0301407;
320
321 /* subtracting admittances from the diagonal elements */
322 bus[9 * i - 1].yff[0] -= 17.3611;
323 bus[9 * i - 1].yff[1] -= -0.0301407;
324
325 bus[9 * i].yff[0] -= 17.3611;
326 bus[9 * i].yff[1] -= -0.0301407;
327 }
328
329 /* read branch data */
330 for (i = 0; i < nc; i++) {
331 for (j = 0; j < NBRANCH; j++) {
332 branch[i * 9 + j].id = i * 9 + j;
333
334 row[0] = edgelist[2 * j] * 2;
335 col[0] = edgelist[2 * j + 1] * 2;
336 col[1] = edgelist[2 * j + 1] * 2 + 1;
337 PetscCall(MatGetValues(Ybus, 1, row, 2, col, branch[i * 9 + j].yft)); /*imaginary part of admittance*/
338 }
339 }
340
341 *pgen = gen;
342 *pexc = exc;
343 *pload = load;
344 *pbus = bus;
345 *pbranch = branch;
346 *pedgelist = edgelist;
347
348 PetscCall(VecRestoreArrayRead(V0, &varr));
349
350 /* Destroy unnecessary stuff */
351 PetscCall(MatDestroy(&Ybus));
352 PetscCall(VecDestroy(&V0));
353 PetscFunctionReturn(PETSC_SUCCESS);
354 }
355
SetInitialGuess(DM networkdm,Vec X)356 PetscErrorCode SetInitialGuess(DM networkdm, Vec X)
357 {
358 Bus *bus;
359 Gen *gen;
360 Exc *exc;
361 PetscInt v, vStart, vEnd, offset;
362 PetscInt key, numComps, j;
363 PetscBool ghostvtex;
364 Vec localX;
365 PetscScalar *xarr;
366 PetscScalar Vr = 0, Vi = 0, Vm = 0, Vm2; /* Terminal voltage variables */
367 PetscScalar IGr, IGi; /* Generator real and imaginary current */
368 PetscScalar Eqp, Edp, delta; /* Generator variables */
369 PetscScalar Efd = 0, RF, VR; /* Exciter variables */
370 PetscScalar Vd, Vq; /* Generator dq axis voltages */
371 PetscScalar Id, Iq; /* Generator dq axis currents */
372 PetscScalar theta; /* Generator phase angle */
373 PetscScalar SE;
374 void *component;
375
376 PetscFunctionBegin;
377 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
378 PetscCall(DMGetLocalVector(networkdm, &localX));
379
380 PetscCall(VecSet(X, 0.0));
381 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
382 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
383
384 PetscCall(VecGetArray(localX, &xarr));
385
386 for (v = vStart; v < vEnd; v++) {
387 PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
388 if (ghostvtex) continue;
389
390 PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));
391 for (j = 0; j < numComps; j++) {
392 PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
393 if (key == 1) {
394 bus = (Bus *)(component);
395
396 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
397 xarr[offset] = bus->vr;
398 xarr[offset + 1] = bus->vi;
399
400 Vr = bus->vr;
401 Vi = bus->vi;
402 } else if (key == 2) {
403 gen = (Gen *)(component);
404 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
405 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
406 Vm2 = Vm * Vm;
407 /* Real part of gen current */
408 IGr = (Vr * gen->PG + Vi * gen->QG) / Vm2;
409 /* Imaginary part of gen current */
410 IGi = (Vi * gen->PG - Vr * gen->QG) / Vm2;
411
412 /* Machine angle */
413 delta = atan2(Vi + gen->Xq * IGr, Vr - gen->Xq * IGi);
414 theta = PETSC_PI / 2.0 - delta;
415
416 /* d-axis stator current */
417 Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta);
418
419 /* q-axis stator current */
420 Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta);
421
422 Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
423 Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
424
425 /* d-axis transient EMF */
426 Edp = Vd + gen->Rs * Id - gen->Xqp * Iq;
427
428 /* q-axis transient EMF */
429 Eqp = Vq + gen->Rs * Iq + gen->Xdp * Id;
430
431 gen->TM = gen->PG;
432
433 xarr[offset] = Eqp;
434 xarr[offset + 1] = Edp;
435 xarr[offset + 2] = delta;
436 xarr[offset + 3] = W_S;
437 xarr[offset + 4] = Id;
438 xarr[offset + 5] = Iq;
439
440 Efd = Eqp + (gen->Xd - gen->Xdp) * Id;
441
442 } else if (key == 3) {
443 exc = (Exc *)(component);
444 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
445
446 SE = exc->k1 * PetscExpScalar(exc->k2 * Efd);
447 VR = exc->KE * Efd + SE;
448 RF = exc->KF * Efd / exc->TF;
449
450 xarr[offset] = Efd;
451 xarr[offset + 1] = RF;
452 xarr[offset + 2] = VR;
453
454 exc->Vref = Vm + (VR / exc->KA);
455 }
456 }
457 }
458 PetscCall(VecRestoreArray(localX, &xarr));
459 PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
460 PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
461 PetscCall(DMRestoreLocalVector(networkdm, &localX));
462 PetscFunctionReturn(PETSC_SUCCESS);
463 }
464
465 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar * Fr,PetscScalar * Fi)466 PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
467 {
468 PetscFunctionBegin;
469 *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
470 *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
471 PetscFunctionReturn(PETSC_SUCCESS);
472 }
473
474 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar * Fd,PetscScalar * Fq)475 PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
476 {
477 PetscFunctionBegin;
478 *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
479 *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
480 PetscFunctionReturn(PETSC_SUCCESS);
481 }
482
483 /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx * user)484 PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
485 {
486 DM networkdm;
487 Vec localX, localXdot, localF;
488 PetscInt vfrom, vto, offsetfrom, offsetto;
489 PetscInt v, vStart, vEnd, e;
490 PetscScalar *farr;
491 PetscScalar Vd = 0, Vq = 0, SE;
492 const PetscScalar *xarr, *xdotarr;
493 void *component;
494 PetscScalar Vr = 0, Vi = 0;
495
496 PetscFunctionBegin;
497 PetscCall(VecSet(F, 0.0));
498
499 PetscCall(TSGetDM(ts, &networkdm));
500 PetscCall(DMGetLocalVector(networkdm, &localF));
501 PetscCall(DMGetLocalVector(networkdm, &localX));
502 PetscCall(DMGetLocalVector(networkdm, &localXdot));
503 PetscCall(VecSet(localF, 0.0));
504
505 /* update ghost values of localX and localXdot */
506 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
507 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
508
509 PetscCall(DMGlobalToLocalBegin(networkdm, Xdot, INSERT_VALUES, localXdot));
510 PetscCall(DMGlobalToLocalEnd(networkdm, Xdot, INSERT_VALUES, localXdot));
511
512 PetscCall(VecGetArrayRead(localX, &xarr));
513 PetscCall(VecGetArrayRead(localXdot, &xdotarr));
514 PetscCall(VecGetArray(localF, &farr));
515
516 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
517
518 for (v = vStart; v < vEnd; v++) {
519 PetscInt i, j, offsetbus, offsetgen, offsetexc, key;
520 Bus *bus;
521 Gen *gen;
522 Exc *exc;
523 Load *load;
524 PetscBool ghostvtex;
525 PetscInt numComps;
526 PetscScalar Yffr, Yffi; /* Real and imaginary fault admittances */
527 PetscScalar Vm, Vm2, Vm0;
528 PetscScalar Vr0 = 0, Vi0 = 0;
529 PetscScalar PD, QD;
530
531 PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
532 PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));
533
534 for (j = 0; j < numComps; j++) {
535 PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
536 if (key == 1) {
537 PetscInt nconnedges;
538 const PetscInt *connedges;
539
540 bus = (Bus *)(component);
541 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus));
542 if (!ghostvtex) {
543 Vr = xarr[offsetbus];
544 Vi = xarr[offsetbus + 1];
545
546 Yffr = bus->yff[1];
547 Yffi = bus->yff[0];
548
549 if (user->alg_flg) {
550 Yffr += user->ybusfault[bus->id * 2 + 1];
551 Yffi += user->ybusfault[bus->id * 2];
552 }
553 Vr0 = bus->vr;
554 Vi0 = bus->vi;
555
556 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
557 The generator current injection, IG, and load current injection, ID are added later
558 */
559 farr[offsetbus] += Yffi * Vr + Yffr * Vi; /* imaginary current due to diagonal elements */
560 farr[offsetbus + 1] += Yffr * Vr - Yffi * Vi; /* real current due to diagonal elements */
561 }
562
563 PetscCall(DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges));
564
565 for (i = 0; i < nconnedges; i++) {
566 Branch *branch;
567 PetscInt keye;
568 PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
569 const PetscInt *cone;
570
571 e = connedges[i];
572 PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));
573
574 Yfti = branch->yft[0];
575 Yftr = branch->yft[1];
576
577 PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
578
579 vfrom = cone[0];
580 vto = cone[1];
581
582 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
583 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));
584
585 /* From bus and to bus real and imaginary voltages */
586 Vfr = xarr[offsetfrom];
587 Vfi = xarr[offsetfrom + 1];
588 Vtr = xarr[offsetto];
589 Vti = xarr[offsetto + 1];
590
591 if (vfrom == v) {
592 farr[offsetfrom] += Yftr * Vti + Yfti * Vtr;
593 farr[offsetfrom + 1] += Yftr * Vtr - Yfti * Vti;
594 } else {
595 farr[offsetto] += Yftr * Vfi + Yfti * Vfr;
596 farr[offsetto + 1] += Yftr * Vfr - Yfti * Vfi;
597 }
598 }
599 } else if (key == 2) {
600 if (!ghostvtex) {
601 PetscScalar Eqp, Edp, delta, w; /* Generator variables */
602 PetscScalar Efd; /* Exciter field voltage */
603 PetscScalar Id, Iq; /* Generator dq axis currents */
604 PetscScalar IGr, IGi, Zdq_inv[4], det;
605 PetscScalar Xd, Xdp, Td0p, Xq, Xqp, Tq0p, TM, D, M, Rs; /* Generator parameters */
606
607 gen = (Gen *)(component);
608 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));
609
610 /* Generator state variables */
611 Eqp = xarr[offsetgen];
612 Edp = xarr[offsetgen + 1];
613 delta = xarr[offsetgen + 2];
614 w = xarr[offsetgen + 3];
615 Id = xarr[offsetgen + 4];
616 Iq = xarr[offsetgen + 5];
617
618 /* Generator parameters */
619 Xd = gen->Xd;
620 Xdp = gen->Xdp;
621 Td0p = gen->Td0p;
622 Xq = gen->Xq;
623 Xqp = gen->Xqp;
624 Tq0p = gen->Tq0p;
625 TM = gen->TM;
626 D = gen->D;
627 M = gen->M;
628 Rs = gen->Rs;
629
630 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, 2, &offsetexc));
631 Efd = xarr[offsetexc];
632
633 /* Generator differential equations */
634 farr[offsetgen] = (Eqp + (Xd - Xdp) * Id - Efd) / Td0p + xdotarr[offsetgen];
635 farr[offsetgen + 1] = (Edp - (Xq - Xqp) * Iq) / Tq0p + xdotarr[offsetgen + 1];
636 farr[offsetgen + 2] = -w + W_S + xdotarr[offsetgen + 2];
637 farr[offsetgen + 3] = (-TM + Edp * Id + Eqp * Iq + (Xqp - Xdp) * Id * Iq + D * (w - W_S)) / M + xdotarr[offsetgen + 3];
638
639 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
640
641 /* Algebraic equations for stator currents */
642 det = Rs * Rs + Xdp * Xqp;
643
644 Zdq_inv[0] = Rs / det;
645 Zdq_inv[1] = Xqp / det;
646 Zdq_inv[2] = -Xdp / det;
647 Zdq_inv[3] = Rs / det;
648
649 farr[offsetgen + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
650 farr[offsetgen + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
651
652 PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
653
654 /* Add generator current injection to network */
655 farr[offsetbus] -= IGi;
656 farr[offsetbus + 1] -= IGr;
657 }
658 } else if (key == 3) {
659 if (!ghostvtex) {
660 PetscScalar k1, k2, KE, TE, TF, KA, KF, Vref, TA; /* Generator parameters */
661 PetscScalar Efd, RF, VR; /* Exciter variables */
662
663 exc = (Exc *)(component);
664 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc));
665
666 Efd = xarr[offsetexc];
667 RF = xarr[offsetexc + 1];
668 VR = xarr[offsetexc + 2];
669
670 k1 = exc->k1;
671 k2 = exc->k2;
672 KE = exc->KE;
673 TE = exc->TE;
674 TF = exc->TF;
675 KA = exc->KA;
676 KF = exc->KF;
677 Vref = exc->Vref;
678 TA = exc->TA;
679
680 Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
681 SE = k1 * PetscExpScalar(k2 * Efd);
682
683 /* Exciter differential equations */
684 farr[offsetexc] = (KE * Efd + SE - VR) / TE + xdotarr[offsetexc];
685 farr[offsetexc + 1] = (RF - KF * Efd / TF) / TF + xdotarr[offsetexc + 1];
686 farr[offsetexc + 2] = (VR - KA * RF + KA * KF * Efd / TF - KA * (Vref - Vm)) / TA + xdotarr[offsetexc + 2];
687 }
688 } else if (key == 4) {
689 if (!ghostvtex) {
690 PetscInt k;
691 PetscInt ld_nsegsp;
692 PetscInt ld_nsegsq;
693 PetscScalar *ld_alphap;
694 PetscScalar *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;
695
696 load = (Load *)(component);
697
698 /* Load Parameters */
699 ld_nsegsp = load->ld_nsegsp;
700 ld_alphap = load->ld_alphap;
701 ld_betap = load->ld_betap;
702 ld_nsegsq = load->ld_nsegsq;
703 ld_alphaq = load->ld_alphaq;
704 ld_betaq = load->ld_betaq;
705 PD0 = load->PD0;
706 QD0 = load->QD0;
707
708 Vr = xarr[offsetbus]; /* Real part of generator terminal voltage */
709 Vi = xarr[offsetbus + 1]; /* Imaginary part of the generator terminal voltage */
710 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
711 Vm2 = Vm * Vm;
712 Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
713 PD = QD = 0.0;
714 for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar(Vm / Vm0, ld_betap[k]);
715 for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
716
717 /* Load currents */
718 IDr = (PD * Vr + QD * Vi) / Vm2;
719 IDi = (-QD * Vr + PD * Vi) / Vm2;
720
721 /* Load current contribution to the network */
722 farr[offsetbus] += IDi;
723 farr[offsetbus + 1] += IDr;
724 }
725 }
726 }
727 }
728
729 PetscCall(VecRestoreArrayRead(localX, &xarr));
730 PetscCall(VecRestoreArrayRead(localXdot, &xdotarr));
731 PetscCall(VecRestoreArray(localF, &farr));
732 PetscCall(DMRestoreLocalVector(networkdm, &localX));
733 PetscCall(DMRestoreLocalVector(networkdm, &localXdot));
734
735 PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
736 PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
737 PetscCall(DMRestoreLocalVector(networkdm, &localF));
738 PetscFunctionReturn(PETSC_SUCCESS);
739 }
740
741 /* This function is used for solving the algebraic system only during fault on and
742 off times. It computes the entire F and then zeros out the part corresponding to
743 differential equations
744 F = [0;g(y)];
745 */
AlgFunction(SNES snes,Vec X,Vec F,PetscCtx ctx)746 PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, PetscCtx ctx)
747 {
748 DM networkdm;
749 Vec localX, localF;
750 PetscInt vfrom, vto, offsetfrom, offsetto;
751 PetscInt v, vStart, vEnd, e;
752 PetscScalar *farr;
753 Userctx *user = (Userctx *)ctx;
754 const PetscScalar *xarr;
755 void *component;
756 PetscScalar Vr = 0, Vi = 0;
757
758 PetscFunctionBegin;
759 PetscCall(VecSet(F, 0.0));
760 PetscCall(SNESGetDM(snes, &networkdm));
761 PetscCall(DMGetLocalVector(networkdm, &localF));
762 PetscCall(DMGetLocalVector(networkdm, &localX));
763 PetscCall(VecSet(localF, 0.0));
764
765 /* update ghost values of locaX and locaXdot */
766 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
767 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
768
769 PetscCall(VecGetArrayRead(localX, &xarr));
770 PetscCall(VecGetArray(localF, &farr));
771
772 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
773
774 for (v = vStart; v < vEnd; v++) {
775 PetscInt i, j, offsetbus, offsetgen, key, numComps;
776 PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0 = 0, Vi0 = 0, PD, QD;
777 Bus *bus;
778 Gen *gen;
779 Load *load;
780 PetscBool ghostvtex;
781
782 PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
783 PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));
784
785 for (j = 0; j < numComps; j++) {
786 PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
787 if (key == 1) {
788 PetscInt nconnedges;
789 const PetscInt *connedges;
790
791 bus = (Bus *)(component);
792 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus));
793 if (!ghostvtex) {
794 Vr = xarr[offsetbus];
795 Vi = xarr[offsetbus + 1];
796
797 Yffr = bus->yff[1];
798 Yffi = bus->yff[0];
799 if (user->alg_flg) {
800 Yffr += user->ybusfault[bus->id * 2 + 1];
801 Yffi += user->ybusfault[bus->id * 2];
802 }
803 Vr0 = bus->vr;
804 Vi0 = bus->vi;
805
806 farr[offsetbus] += Yffi * Vr + Yffr * Vi;
807 farr[offsetbus + 1] += Yffr * Vr - Yffi * Vi;
808 }
809 PetscCall(DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges));
810
811 for (i = 0; i < nconnedges; i++) {
812 Branch *branch;
813 PetscInt keye;
814 PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
815 const PetscInt *cone;
816
817 e = connedges[i];
818 PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));
819
820 Yfti = branch->yft[0];
821 Yftr = branch->yft[1];
822
823 PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
824 vfrom = cone[0];
825 vto = cone[1];
826
827 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
828 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));
829
830 /*From bus and to bus real and imaginary voltages */
831 Vfr = xarr[offsetfrom];
832 Vfi = xarr[offsetfrom + 1];
833 Vtr = xarr[offsetto];
834 Vti = xarr[offsetto + 1];
835
836 if (vfrom == v) {
837 farr[offsetfrom] += Yftr * Vti + Yfti * Vtr;
838 farr[offsetfrom + 1] += Yftr * Vtr - Yfti * Vti;
839 } else {
840 farr[offsetto] += Yftr * Vfi + Yfti * Vfr;
841 farr[offsetto + 1] += Yftr * Vfr - Yfti * Vfi;
842 }
843 }
844 } else if (key == 2) {
845 if (!ghostvtex) {
846 PetscScalar Eqp, Edp, delta; /* Generator variables */
847 PetscScalar Id, Iq; /* Generator dq axis currents */
848 PetscScalar Vd, Vq, IGr, IGi, Zdq_inv[4], det;
849 PetscScalar Xdp, Xqp, Rs; /* Generator parameters */
850
851 gen = (Gen *)(component);
852 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));
853
854 /* Generator state variables */
855 Eqp = xarr[offsetgen];
856 Edp = xarr[offsetgen + 1];
857 delta = xarr[offsetgen + 2];
858 /* w = xarr[idx+3]; not being used */
859 Id = xarr[offsetgen + 4];
860 Iq = xarr[offsetgen + 5];
861
862 /* Generator parameters */
863 Xdp = gen->Xdp;
864 Xqp = gen->Xqp;
865 Rs = gen->Rs;
866
867 /* Set generator differential equation residual functions to zero */
868 farr[offsetgen] = 0;
869 farr[offsetgen + 1] = 0;
870 farr[offsetgen + 2] = 0;
871 farr[offsetgen + 3] = 0;
872
873 PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
874
875 /* Algebraic equations for stator currents */
876 det = Rs * Rs + Xdp * Xqp;
877
878 Zdq_inv[0] = Rs / det;
879 Zdq_inv[1] = Xqp / det;
880 Zdq_inv[2] = -Xdp / det;
881 Zdq_inv[3] = Rs / det;
882
883 farr[offsetgen + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
884 farr[offsetgen + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
885
886 /* Add generator current injection to network */
887 PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
888
889 farr[offsetbus] -= IGi;
890 farr[offsetbus + 1] -= IGr;
891
892 /* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */
893 }
894 } else if (key == 3) {
895 if (!ghostvtex) {
896 PetscInt offsetexc;
897 PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc));
898 /* Set exciter differential equation residual functions equal to zero*/
899 farr[offsetexc] = 0;
900 farr[offsetexc + 1] = 0;
901 farr[offsetexc + 2] = 0;
902 }
903 } else if (key == 4) {
904 if (!ghostvtex) {
905 PetscInt k, ld_nsegsp, ld_nsegsq;
906 PetscScalar *ld_alphap, *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;
907
908 load = (Load *)(component);
909
910 /* Load Parameters */
911 ld_nsegsp = load->ld_nsegsp;
912 ld_alphap = load->ld_alphap;
913 ld_betap = load->ld_betap;
914 ld_nsegsq = load->ld_nsegsq;
915 ld_alphaq = load->ld_alphaq;
916 ld_betaq = load->ld_betaq;
917
918 PD0 = load->PD0;
919 QD0 = load->QD0;
920
921 Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
922 Vm2 = Vm * Vm;
923 Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
924 PD = QD = 0.0;
925 for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar(Vm / Vm0, ld_betap[k]);
926 for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
927
928 /* Load currents */
929 IDr = (PD * Vr + QD * Vi) / Vm2;
930 IDi = (-QD * Vr + PD * Vi) / Vm2;
931
932 farr[offsetbus] += IDi;
933 farr[offsetbus + 1] += IDr;
934 }
935 }
936 }
937 }
938
939 PetscCall(VecRestoreArrayRead(localX, &xarr));
940 PetscCall(VecRestoreArray(localF, &farr));
941 PetscCall(DMRestoreLocalVector(networkdm, &localX));
942
943 PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
944 PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
945 PetscCall(DMRestoreLocalVector(networkdm, &localF));
946 PetscFunctionReturn(PETSC_SUCCESS);
947 }
948
main(int argc,char ** argv)949 int main(int argc, char **argv)
950 {
951 PetscInt i, j, *edgelist = NULL, eStart, eEnd, vStart, vEnd;
952 PetscInt genj, excj, loadj, componentkey[5];
953 PetscInt nc = 1; /* No. of copies (default = 1) */
954 PetscMPIInt size, rank;
955 Vec X, F_alg;
956 TS ts;
957 SNES snes_alg, snes;
958 Bus *bus;
959 Branch *branch;
960 Gen *gen;
961 Exc *exc;
962 Load *load;
963 DM networkdm;
964 PetscLogStage stage1;
965 Userctx user;
966 KSP ksp;
967 PC pc;
968 PetscInt numEdges = 0;
969
970 PetscFunctionBeginUser;
971 PetscCall(PetscInitialize(&argc, &argv, "ex9busnetworkops", help));
972 PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
973 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
974 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
975
976 /* Read initial voltage vector and Ybus */
977 if (rank == 0) PetscCall(read_data(nc, &gen, &exc, &load, &bus, &branch, &edgelist));
978
979 PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
980 PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(Branch), &componentkey[0]));
981 PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(Bus), &componentkey[1]));
982 PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(Gen), &componentkey[2]));
983 PetscCall(DMNetworkRegisterComponent(networkdm, "excstruct", sizeof(Exc), &componentkey[3]));
984 PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(Load), &componentkey[4]));
985
986 PetscCall(PetscLogStageRegister("Create network", &stage1));
987 PetscCall(PetscLogStagePush(stage1));
988
989 /* Set local number of edges and edge connectivity */
990 if (rank == 0) numEdges = NBRANCH * nc + (nc - 1);
991 PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
992 PetscCall(DMNetworkAddSubnetwork(networkdm, NULL, numEdges, edgelist, NULL));
993
994 /* Set up the network layout */
995 PetscCall(DMNetworkLayoutSetUp(networkdm));
996
997 if (rank == 0) PetscCall(PetscFree(edgelist));
998
999 /* Add network components (physical parameters of nodes and branches) and number of variables */
1000 if (rank == 0) {
1001 PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
1002 genj = 0;
1003 loadj = 0;
1004 excj = 0;
1005 for (i = eStart; i < eEnd; i++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[0], &branch[i - eStart], 0));
1006
1007 PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
1008
1009 for (i = vStart; i < vEnd; i++) {
1010 PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[1], &bus[i - vStart], 2));
1011 if (bus[i - vStart].nofgen) {
1012 for (j = 0; j < bus[i - vStart].nofgen; j++) {
1013 /* Add generator */
1014 PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[2], &gen[genj++], 6));
1015 /* Add exciter */
1016 PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[3], &exc[excj++], 3));
1017 }
1018 }
1019 if (bus[i - vStart].nofload) {
1020 for (j = 0; j < bus[i - vStart].nofload; j++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[4], &load[loadj++], 0));
1021 }
1022 }
1023 }
1024
1025 PetscCall(DMSetUp(networkdm));
1026
1027 if (rank == 0) PetscCall(PetscFree5(bus, gen, load, branch, exc));
1028
1029 /* for parallel options: Network partitioning and distribution of data */
1030 if (size > 1) PetscCall(DMNetworkDistribute(&networkdm, 0));
1031 PetscCall(PetscLogStagePop());
1032
1033 PetscCall(DMCreateGlobalVector(networkdm, &X));
1034
1035 PetscCall(SetInitialGuess(networkdm, X));
1036
1037 /* Options for fault simulation */
1038 PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1039 user.tfaulton = 0.02;
1040 user.tfaultoff = 0.05;
1041 user.Rfault = 0.0001;
1042 user.faultbus = 8;
1043 PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
1044 PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
1045 PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
1046 user.t0 = 0.0;
1047 user.tmax = 0.1;
1048 PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
1049 PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
1050
1051 PetscCall(PetscMalloc1(18 * nc, &user.ybusfault));
1052 for (i = 0; i < 18 * nc; i++) user.ybusfault[i] = 0;
1053 user.ybusfault[user.faultbus * 2 + 1] = 1 / user.Rfault;
1054 PetscOptionsEnd();
1055
1056 /* Setup TS solver */
1057 /*--------------------------------------------------------*/
1058 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1059 PetscCall(TSSetDM(ts, (DM)networkdm));
1060 PetscCall(TSSetType(ts, TSCN));
1061
1062 PetscCall(TSGetSNES(ts, &snes));
1063 PetscCall(SNESGetKSP(snes, &ksp));
1064 PetscCall(KSPGetPC(ksp, &pc));
1065 PetscCall(PCSetType(pc, PCBJACOBI));
1066
1067 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));
1068 PetscCall(TSSetMaxTime(ts, user.tfaulton));
1069 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1070 PetscCall(TSSetTimeStep(ts, 0.01));
1071 PetscCall(TSSetFromOptions(ts));
1072
1073 /*user.alg_flg = PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix.
1074 eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */
1075 user.alg_flg = PETSC_FALSE;
1076
1077 /* Prefault period */
1078 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "... (1) Prefault period ... \n"));
1079
1080 PetscCall(TSSetSolution(ts, X));
1081 PetscCall(TSSetUp(ts));
1082 PetscCall(TSSolve(ts, X));
1083
1084 /* Create the nonlinear solver for solving the algebraic system */
1085 PetscCall(VecDuplicate(X, &F_alg));
1086
1087 PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1088 PetscCall(SNESSetDM(snes_alg, (DM)networkdm));
1089 PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));
1090 PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1091 PetscCall(SNESSetFromOptions(snes_alg));
1092
1093 /* Apply disturbance - resistive fault at user.faultbus */
1094 /* This is done by adding shunt conductance to the diagonal location
1095 in the Ybus matrix */
1096 user.alg_flg = PETSC_TRUE;
1097
1098 /* Solve the algebraic equations */
1099 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (2) Apply disturbance, solve algebraic equations ... \n"));
1100 PetscCall(SNESSolve(snes_alg, NULL, X));
1101
1102 /* Disturbance period */
1103 PetscCall(TSSetTime(ts, user.tfaulton));
1104 PetscCall(TSSetMaxTime(ts, user.tfaultoff));
1105 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1106 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));
1107
1108 user.alg_flg = PETSC_TRUE;
1109 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (3) Disturbance period ... \n"));
1110 PetscCall(TSSolve(ts, X));
1111
1112 /* Remove the fault */
1113 PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));
1114
1115 user.alg_flg = PETSC_FALSE;
1116 /* Solve the algebraic equations */
1117 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (4) Remove fault, solve algebraic equations ... \n"));
1118 PetscCall(SNESSolve(snes_alg, NULL, X));
1119 PetscCall(SNESDestroy(&snes_alg));
1120
1121 /* Post-disturbance period */
1122 PetscCall(TSSetTime(ts, user.tfaultoff));
1123 PetscCall(TSSetMaxTime(ts, user.tmax));
1124 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1125 PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));
1126
1127 user.alg_flg = PETSC_FALSE;
1128 if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (5) Post-disturbance period ... \n"));
1129 PetscCall(TSSolve(ts, X));
1130
1131 PetscCall(PetscFree(user.ybusfault));
1132 PetscCall(VecDestroy(&F_alg));
1133 PetscCall(VecDestroy(&X));
1134 PetscCall(DMDestroy(&networkdm));
1135 PetscCall(TSDestroy(&ts));
1136 PetscCall(PetscFinalize());
1137 return 0;
1138 }
1139
1140 /*TEST
1141
1142 build:
1143 requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1144
1145 test:
1146 args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1147 localrunfiles: X.bin Ybus.bin ex9busnetworkops
1148
1149 test:
1150 suffix: 2
1151 nsize: 2
1152 args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1153 localrunfiles: X.bin Ybus.bin ex9busnetworkops
1154
1155 TEST*/
1156