Actual source code: ex9busdmnetwork.c
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";
9: /*
10: This example was modified from ex9busdmnetwork.c.
11: */
13: #include <petscts.h>
14: #include <petscdmnetwork.h>
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 */
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 */
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;
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;
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;
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
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;
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;
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;
103: /* Used to read data into the DMNetwork components */
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;
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};
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 */
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) */
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;
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));
164: PetscCall(MatCreate(PETSC_COMM_SELF, &Ybus));
165: PetscCall(MatSetType(Ybus, MATBAIJ));
166: PetscCall(MatLoad(Ybus, Ybusview));
168: /* Destroy unnecessary stuff */
169: PetscCall(PetscViewerDestroy(&Xview));
170: PetscCall(PetscViewerDestroy(&Ybusview));
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");
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];
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));
186: PetscCall(VecGetArrayRead(V0, &varr));
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: }
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: }
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: }
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];
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];
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];
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];
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));
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: }
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;
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;
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;
325: bus[9 * i].yff[0] -= 17.3611;
326: bus[9 * i].yff[1] -= -0.0301407;
327: }
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;
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: }
341: *pgen = gen;
342: *pexc = exc;
343: *pload = load;
344: *pbus = bus;
345: *pbranch = branch;
346: *pedgelist = edgelist;
348: PetscCall(VecRestoreArrayRead(V0, &varr));
350: /* Destroy unnecessary stuff */
351: PetscCall(MatDestroy(&Ybus));
352: PetscCall(VecDestroy(&V0));
353: PetscFunctionReturn(PETSC_SUCCESS);
354: }
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;
376: PetscFunctionBegin;
377: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
378: PetscCall(DMGetLocalVector(networkdm, &localX));
380: PetscCall(VecSet(X, 0.0));
381: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
382: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
384: PetscCall(VecGetArray(localX, &xarr));
386: for (v = vStart; v < vEnd; v++) {
387: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
388: if (ghostvtex) continue;
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);
396: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
397: xarr[offset] = bus->vr;
398: xarr[offset + 1] = bus->vi;
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;
412: /* Machine angle */
413: delta = atan2(Vi + gen->Xq * IGr, Vr - gen->Xq * IGi);
414: theta = PETSC_PI / 2.0 - delta;
416: /* d-axis stator current */
417: Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta);
419: /* q-axis stator current */
420: Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta);
422: Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
423: Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
425: /* d-axis transient EMF */
426: Edp = Vd + gen->Rs * Id - gen->Xqp * Iq;
428: /* q-axis transient EMF */
429: Eqp = Vq + gen->Rs * Iq + gen->Xdp * Id;
431: gen->TM = gen->PG;
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;
440: Efd = Eqp + (gen->Xd - gen->Xdp) * Id;
442: } else if (key == 3) {
443: exc = (Exc *)(component);
444: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
446: SE = exc->k1 * PetscExpScalar(exc->k2 * Efd);
447: VR = exc->KE * Efd + SE;
448: RF = exc->KF * Efd / exc->TF;
450: xarr[offset] = Efd;
451: xarr[offset + 1] = RF;
452: xarr[offset + 2] = VR;
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: }
465: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
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: }
474: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
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: }
483: /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
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;
496: PetscFunctionBegin;
497: PetscCall(VecSet(F, 0.0));
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));
505: /* update ghost values of localX and localXdot */
506: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
507: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
509: PetscCall(DMGlobalToLocalBegin(networkdm, Xdot, INSERT_VALUES, localXdot));
510: PetscCall(DMGlobalToLocalEnd(networkdm, Xdot, INSERT_VALUES, localXdot));
512: PetscCall(VecGetArrayRead(localX, &xarr));
513: PetscCall(VecGetArrayRead(localXdot, &xdotarr));
514: PetscCall(VecGetArray(localF, &farr));
516: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
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;
531: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
532: PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));
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;
540: bus = (Bus *)(component);
541: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus));
542: if (!ghostvtex) {
543: Vr = xarr[offsetbus];
544: Vi = xarr[offsetbus + 1];
546: Yffr = bus->yff[1];
547: Yffi = bus->yff[0];
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;
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: }
563: PetscCall(DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges));
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;
571: e = connedges[i];
572: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));
574: Yfti = branch->yft[0];
575: Yftr = branch->yft[1];
577: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
579: vfrom = cone[0];
580: vto = cone[1];
582: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
583: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));
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];
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 */
607: gen = (Gen *)(component);
608: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));
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];
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;
630: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, 2, &offsetexc));
631: Efd = xarr[offsetexc];
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];
639: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
641: /* Algebraic equations for stator currents */
642: det = Rs * Rs + Xdp * Xqp;
644: Zdq_inv[0] = Rs / det;
645: Zdq_inv[1] = Xqp / det;
646: Zdq_inv[2] = -Xdp / det;
647: Zdq_inv[3] = Rs / det;
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;
652: PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
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 */
663: exc = (Exc *)(component);
664: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc));
666: Efd = xarr[offsetexc];
667: RF = xarr[offsetexc + 1];
668: VR = xarr[offsetexc + 2];
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;
680: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
681: SE = k1 * PetscExpScalar(k2 * Efd);
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 ld_nsegsp;
691: PetscInt ld_nsegsq;
692: PetscScalar *ld_alphap;
693: PetscScalar *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;
695: load = (Load *)(component);
697: /* Load Parameters */
698: ld_nsegsp = load->ld_nsegsp;
699: ld_alphap = load->ld_alphap;
700: ld_betap = load->ld_betap;
701: ld_nsegsq = load->ld_nsegsq;
702: ld_alphaq = load->ld_alphaq;
703: ld_betaq = load->ld_betaq;
704: PD0 = load->PD0;
705: QD0 = load->QD0;
707: Vr = xarr[offsetbus]; /* Real part of generator terminal voltage */
708: Vi = xarr[offsetbus + 1]; /* Imaginary part of the generator terminal voltage */
709: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
710: Vm2 = Vm * Vm;
711: Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
712: PD = QD = 0.0;
713: for (PetscInt k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar(Vm / Vm0, ld_betap[k]);
714: for (PetscInt k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
716: /* Load currents */
717: IDr = (PD * Vr + QD * Vi) / Vm2;
718: IDi = (-QD * Vr + PD * Vi) / Vm2;
720: /* Load current contribution to the network */
721: farr[offsetbus] += IDi;
722: farr[offsetbus + 1] += IDr;
723: }
724: }
725: }
726: }
728: PetscCall(VecRestoreArrayRead(localX, &xarr));
729: PetscCall(VecRestoreArrayRead(localXdot, &xdotarr));
730: PetscCall(VecRestoreArray(localF, &farr));
731: PetscCall(DMRestoreLocalVector(networkdm, &localX));
732: PetscCall(DMRestoreLocalVector(networkdm, &localXdot));
734: PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
735: PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
736: PetscCall(DMRestoreLocalVector(networkdm, &localF));
737: PetscFunctionReturn(PETSC_SUCCESS);
738: }
740: /* This function is used for solving the algebraic system only during fault on and
741: off times. It computes the entire F and then zeros out the part corresponding to
742: differential equations
743: F = [0;g(y)];
744: */
745: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, PetscCtx ctx)
746: {
747: DM networkdm;
748: Vec localX, localF;
749: PetscInt vfrom, vto, offsetfrom, offsetto;
750: PetscInt v, vStart, vEnd, e;
751: PetscScalar *farr;
752: Userctx *user = (Userctx *)ctx;
753: const PetscScalar *xarr;
754: void *component;
755: PetscScalar Vr = 0, Vi = 0;
757: PetscFunctionBegin;
758: PetscCall(VecSet(F, 0.0));
759: PetscCall(SNESGetDM(snes, &networkdm));
760: PetscCall(DMGetLocalVector(networkdm, &localF));
761: PetscCall(DMGetLocalVector(networkdm, &localX));
762: PetscCall(VecSet(localF, 0.0));
764: /* update ghost values of locaX and locaXdot */
765: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
766: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
768: PetscCall(VecGetArrayRead(localX, &xarr));
769: PetscCall(VecGetArray(localF, &farr));
771: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
773: for (v = vStart; v < vEnd; v++) {
774: PetscInt i, j, offsetbus, offsetgen, key, numComps;
775: PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0 = 0, Vi0 = 0, PD, QD;
776: Bus *bus;
777: Gen *gen;
778: Load *load;
779: PetscBool ghostvtex;
781: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
782: PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));
784: for (j = 0; j < numComps; j++) {
785: PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
786: if (key == 1) {
787: PetscInt nconnedges;
788: const PetscInt *connedges;
790: bus = (Bus *)(component);
791: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus));
792: if (!ghostvtex) {
793: Vr = xarr[offsetbus];
794: Vi = xarr[offsetbus + 1];
796: Yffr = bus->yff[1];
797: Yffi = bus->yff[0];
798: if (user->alg_flg) {
799: Yffr += user->ybusfault[bus->id * 2 + 1];
800: Yffi += user->ybusfault[bus->id * 2];
801: }
802: Vr0 = bus->vr;
803: Vi0 = bus->vi;
805: farr[offsetbus] += Yffi * Vr + Yffr * Vi;
806: farr[offsetbus + 1] += Yffr * Vr - Yffi * Vi;
807: }
808: PetscCall(DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges));
810: for (i = 0; i < nconnedges; i++) {
811: Branch *branch;
812: PetscInt keye;
813: PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
814: const PetscInt *cone;
816: e = connedges[i];
817: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));
819: Yfti = branch->yft[0];
820: Yftr = branch->yft[1];
822: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
823: vfrom = cone[0];
824: vto = cone[1];
826: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
827: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));
829: /*From bus and to bus real and imaginary voltages */
830: Vfr = xarr[offsetfrom];
831: Vfi = xarr[offsetfrom + 1];
832: Vtr = xarr[offsetto];
833: Vti = xarr[offsetto + 1];
835: if (vfrom == v) {
836: farr[offsetfrom] += Yftr * Vti + Yfti * Vtr;
837: farr[offsetfrom + 1] += Yftr * Vtr - Yfti * Vti;
838: } else {
839: farr[offsetto] += Yftr * Vfi + Yfti * Vfr;
840: farr[offsetto + 1] += Yftr * Vfr - Yfti * Vfi;
841: }
842: }
843: } else if (key == 2) {
844: if (!ghostvtex) {
845: PetscScalar Eqp, Edp, delta; /* Generator variables */
846: PetscScalar Id, Iq; /* Generator dq axis currents */
847: PetscScalar Vd, Vq, IGr, IGi, Zdq_inv[4], det;
848: PetscScalar Xdp, Xqp, Rs; /* Generator parameters */
850: gen = (Gen *)(component);
851: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));
853: /* Generator state variables */
854: Eqp = xarr[offsetgen];
855: Edp = xarr[offsetgen + 1];
856: delta = xarr[offsetgen + 2];
857: /* w = xarr[idx+3]; not being used */
858: Id = xarr[offsetgen + 4];
859: Iq = xarr[offsetgen + 5];
861: /* Generator parameters */
862: Xdp = gen->Xdp;
863: Xqp = gen->Xqp;
864: Rs = gen->Rs;
866: /* Set generator differential equation residual functions to zero */
867: farr[offsetgen] = 0;
868: farr[offsetgen + 1] = 0;
869: farr[offsetgen + 2] = 0;
870: farr[offsetgen + 3] = 0;
872: PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));
874: /* Algebraic equations for stator currents */
875: det = Rs * Rs + Xdp * Xqp;
877: Zdq_inv[0] = Rs / det;
878: Zdq_inv[1] = Xqp / det;
879: Zdq_inv[2] = -Xdp / det;
880: Zdq_inv[3] = Rs / det;
882: farr[offsetgen + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
883: farr[offsetgen + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
885: /* Add generator current injection to network */
886: PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));
888: farr[offsetbus] -= IGi;
889: farr[offsetbus + 1] -= IGr;
891: /* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */
892: }
893: } else if (key == 3) {
894: if (!ghostvtex) {
895: PetscInt offsetexc;
896: PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc));
897: /* Set exciter differential equation residual functions equal to zero*/
898: farr[offsetexc] = 0;
899: farr[offsetexc + 1] = 0;
900: farr[offsetexc + 2] = 0;
901: }
902: } else if (key == 4) {
903: if (!ghostvtex) {
904: PetscInt k, ld_nsegsp, ld_nsegsq;
905: PetscScalar *ld_alphap, *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;
907: load = (Load *)(component);
909: /* Load Parameters */
910: ld_nsegsp = load->ld_nsegsp;
911: ld_alphap = load->ld_alphap;
912: ld_betap = load->ld_betap;
913: ld_nsegsq = load->ld_nsegsq;
914: ld_alphaq = load->ld_alphaq;
915: ld_betaq = load->ld_betaq;
917: PD0 = load->PD0;
918: QD0 = load->QD0;
920: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
921: Vm2 = Vm * Vm;
922: Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
923: PD = QD = 0.0;
924: for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar(Vm / Vm0, ld_betap[k]);
925: for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar(Vm / Vm0, ld_betaq[k]);
927: /* Load currents */
928: IDr = (PD * Vr + QD * Vi) / Vm2;
929: IDi = (-QD * Vr + PD * Vi) / Vm2;
931: farr[offsetbus] += IDi;
932: farr[offsetbus + 1] += IDr;
933: }
934: }
935: }
936: }
938: PetscCall(VecRestoreArrayRead(localX, &xarr));
939: PetscCall(VecRestoreArray(localF, &farr));
940: PetscCall(DMRestoreLocalVector(networkdm, &localX));
942: PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
943: PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
944: PetscCall(DMRestoreLocalVector(networkdm, &localF));
945: PetscFunctionReturn(PETSC_SUCCESS);
946: }
948: int main(int argc, char **argv)
949: {
950: PetscInt i, j, *edgelist = NULL, eStart, eEnd, vStart, vEnd;
951: PetscInt genj, excj, loadj, componentkey[5];
952: PetscInt nc = 1; /* No. of copies (default = 1) */
953: PetscMPIInt size, rank;
954: Vec X, F_alg;
955: TS ts;
956: SNES snes_alg, snes;
957: Bus *bus;
958: Branch *branch;
959: Gen *gen;
960: Exc *exc;
961: Load *load;
962: DM networkdm;
963: PetscLogStage stage1;
964: Userctx user;
965: KSP ksp;
966: PC pc;
967: PetscInt numEdges = 0;
969: PetscFunctionBeginUser;
970: PetscCall(PetscInitialize(&argc, &argv, "ex9busnetworkops", help));
971: PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
972: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
973: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
975: /* Read initial voltage vector and Ybus */
976: if (rank == 0) PetscCall(read_data(nc, &gen, &exc, &load, &bus, &branch, &edgelist));
978: PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
979: PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(Branch), &componentkey[0]));
980: PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(Bus), &componentkey[1]));
981: PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(Gen), &componentkey[2]));
982: PetscCall(DMNetworkRegisterComponent(networkdm, "excstruct", sizeof(Exc), &componentkey[3]));
983: PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(Load), &componentkey[4]));
985: PetscCall(PetscLogStageRegister("Create network", &stage1));
986: PetscCall(PetscLogStagePush(stage1));
988: /* Set local number of edges and edge connectivity */
989: if (rank == 0) numEdges = NBRANCH * nc + (nc - 1);
990: PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
991: PetscCall(DMNetworkAddSubnetwork(networkdm, NULL, numEdges, edgelist, NULL));
993: /* Set up the network layout */
994: PetscCall(DMNetworkLayoutSetUp(networkdm));
996: if (rank == 0) PetscCall(PetscFree(edgelist));
998: /* Add network components (physical parameters of nodes and branches) and number of variables */
999: if (rank == 0) {
1000: PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
1001: genj = 0;
1002: loadj = 0;
1003: excj = 0;
1004: for (i = eStart; i < eEnd; i++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[0], &branch[i - eStart], 0));
1006: PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
1008: for (i = vStart; i < vEnd; i++) {
1009: PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[1], &bus[i - vStart], 2));
1010: if (bus[i - vStart].nofgen) {
1011: for (j = 0; j < bus[i - vStart].nofgen; j++) {
1012: /* Add generator */
1013: PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[2], &gen[genj++], 6));
1014: /* Add exciter */
1015: PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[3], &exc[excj++], 3));
1016: }
1017: }
1018: if (bus[i - vStart].nofload) {
1019: for (j = 0; j < bus[i - vStart].nofload; j++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[4], &load[loadj++], 0));
1020: }
1021: }
1022: }
1024: PetscCall(DMSetUp(networkdm));
1026: if (rank == 0) PetscCall(PetscFree5(bus, gen, load, branch, exc));
1028: /* for parallel options: Network partitioning and distribution of data */
1029: if (size > 1) PetscCall(DMNetworkDistribute(&networkdm, 0));
1030: PetscCall(PetscLogStagePop());
1032: PetscCall(DMCreateGlobalVector(networkdm, &X));
1034: PetscCall(SetInitialGuess(networkdm, X));
1036: /* Options for fault simulation */
1037: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1038: user.tfaulton = 0.02;
1039: user.tfaultoff = 0.05;
1040: user.Rfault = 0.0001;
1041: user.faultbus = 8;
1042: PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
1043: PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
1044: PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
1045: user.t0 = 0.0;
1046: user.tmax = 0.1;
1047: PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
1048: PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));
1050: PetscCall(PetscMalloc1(18 * nc, &user.ybusfault));
1051: for (i = 0; i < 18 * nc; i++) user.ybusfault[i] = 0;
1052: user.ybusfault[user.faultbus * 2 + 1] = 1 / user.Rfault;
1053: PetscOptionsEnd();
1055: /* Setup TS solver */
1056: /*--------------------------------------------------------*/
1057: PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1058: PetscCall(TSSetDM(ts, (DM)networkdm));
1059: PetscCall(TSSetType(ts, TSCN));
1061: PetscCall(TSGetSNES(ts, &snes));
1062: PetscCall(SNESGetKSP(snes, &ksp));
1063: PetscCall(KSPGetPC(ksp, &pc));
1064: PetscCall(PCSetType(pc, PCBJACOBI));
1066: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));
1067: PetscCall(TSSetMaxTime(ts, user.tfaulton));
1068: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1069: PetscCall(TSSetTimeStep(ts, 0.01));
1070: PetscCall(TSSetFromOptions(ts));
1072: /*user.alg_flg = PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix.
1073: eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */
1074: user.alg_flg = PETSC_FALSE;
1076: /* Prefault period */
1077: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "... (1) Prefault period ... \n"));
1079: PetscCall(TSSetSolution(ts, X));
1080: PetscCall(TSSetUp(ts));
1081: PetscCall(TSSolve(ts, X));
1083: /* Create the nonlinear solver for solving the algebraic system */
1084: PetscCall(VecDuplicate(X, &F_alg));
1086: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1087: PetscCall(SNESSetDM(snes_alg, (DM)networkdm));
1088: PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));
1089: PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1090: PetscCall(SNESSetFromOptions(snes_alg));
1092: /* Apply disturbance - resistive fault at user.faultbus */
1093: /* This is done by adding shunt conductance to the diagonal location
1094: in the Ybus matrix */
1095: user.alg_flg = PETSC_TRUE;
1097: /* Solve the algebraic equations */
1098: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (2) Apply disturbance, solve algebraic equations ... \n"));
1099: PetscCall(SNESSolve(snes_alg, NULL, X));
1101: /* Disturbance period */
1102: PetscCall(TSSetTime(ts, user.tfaulton));
1103: PetscCall(TSSetMaxTime(ts, user.tfaultoff));
1104: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1105: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));
1107: user.alg_flg = PETSC_TRUE;
1108: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (3) Disturbance period ... \n"));
1109: PetscCall(TSSolve(ts, X));
1111: /* Remove the fault */
1112: PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));
1114: user.alg_flg = PETSC_FALSE;
1115: /* Solve the algebraic equations */
1116: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (4) Remove fault, solve algebraic equations ... \n"));
1117: PetscCall(SNESSolve(snes_alg, NULL, X));
1118: PetscCall(SNESDestroy(&snes_alg));
1120: /* Post-disturbance period */
1121: PetscCall(TSSetTime(ts, user.tfaultoff));
1122: PetscCall(TSSetMaxTime(ts, user.tmax));
1123: PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1124: PetscCall(TSSetIFunction(ts, NULL, (TSIFunctionFn *)FormIFunction, &user));
1126: user.alg_flg = PETSC_FALSE;
1127: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (5) Post-disturbance period ... \n"));
1128: PetscCall(TSSolve(ts, X));
1130: PetscCall(PetscFree(user.ybusfault));
1131: PetscCall(VecDestroy(&F_alg));
1132: PetscCall(VecDestroy(&X));
1133: PetscCall(DMDestroy(&networkdm));
1134: PetscCall(TSDestroy(&ts));
1135: PetscCall(PetscFinalize());
1136: return 0;
1137: }
1139: /*TEST
1141: build:
1142: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1144: test:
1145: args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1146: localrunfiles: X.bin Ybus.bin ex9busnetworkops
1148: test:
1149: suffix: 2
1150: nsize: 2
1151: args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1152: localrunfiles: X.bin Ybus.bin ex9busnetworkops
1154: TEST*/