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