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: /* Alocate 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;
655: }
656: } else if (key == 3) {
657: if (!ghostvtex) {
658: PetscScalar k1,k2,KE,TE,TF,KA,KF,Vref,TA; /* Generator parameters */
659: PetscScalar Efd,RF,VR; /* Exciter variables */
661: exc = (Exc*)(component);
662: DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetexc);
664: Efd = xarr[offsetexc];
665: RF = xarr[offsetexc+1];
666: VR = xarr[offsetexc+2];
668: k1 = exc->k1;
669: k2 = exc->k2;
670: KE = exc->KE;
671: TE = exc->TE;
672: TF = exc->TF;
673: KA = exc->KA;
674: KF = exc->KF;
675: Vref = exc->Vref;
676: TA = exc->TA;
678: Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);
679: SE = k1*PetscExpScalar(k2*Efd);
681: /* Exciter differential equations */
682: farr[offsetexc] = (KE*Efd + SE - VR)/TE + xdotarr[offsetexc];
683: farr[offsetexc+1] = (RF - KF*Efd/TF)/TF + xdotarr[offsetexc+1];
684: farr[offsetexc+2] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[offsetexc+2];
686: }
687: } else if (key ==4) {
688: if (!ghostvtex) {
689: PetscInt k;
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 (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]);
714: for (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: VecRestoreArrayRead(localX,&xarr);
729: VecRestoreArrayRead(localXdot,&xdotarr);
730: VecRestoreArray(localF,&farr);
731: DMRestoreLocalVector(networkdm,&localX);
732: DMRestoreLocalVector(networkdm,&localXdot);
734: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
735: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
736: DMRestoreLocalVector(networkdm,&localF);
737: return 0;
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, void *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: VecSet(F,0.0);
758: SNESGetDM(snes,&networkdm);
759: DMGetLocalVector(networkdm,&localF);
760: DMGetLocalVector(networkdm,&localX);
761: VecSet(localF,0.0);
763: /* update ghost values of locaX and locaXdot */
764: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
765: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
767: VecGetArrayRead(localX,&xarr);
768: VecGetArray(localF,&farr);
770: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
772: for (v=vStart; v < vEnd; v++) {
773: PetscInt i,j,offsetbus,offsetgen,key,numComps;
774: PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0=0, Vi0=0, PD, QD;
775: Bus *bus;
776: Gen *gen;
777: Load *load;
778: PetscBool ghostvtex;
780: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
781: DMNetworkGetNumComponents(networkdm,v,&numComps);
783: for (j = 0; j < numComps; j++) {
784: DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL);
785: if (key == 1) {
786: PetscInt nconnedges;
787: const PetscInt *connedges;
789: bus = (Bus*)(component);
790: DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetbus);
791: if (!ghostvtex) {
792: Vr = xarr[offsetbus];
793: Vi = xarr[offsetbus+1];
795: Yffr = bus->yff[1];
796: Yffi = bus->yff[0];
797: if (user->alg_flg) {
798: Yffr += user->ybusfault[bus->id*2+1];
799: Yffi += user->ybusfault[bus->id*2];
800: }
801: Vr0 = bus->vr;
802: Vi0 = bus->vi;
804: farr[offsetbus] += Yffi*Vr + Yffr*Vi;
805: farr[offsetbus+1] += Yffr*Vr - Yffi*Vi;
806: }
807: DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
809: for (i=0; i < nconnedges; i++) {
810: Branch *branch;
811: PetscInt keye;
812: PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
813: const PetscInt *cone;
815: e = connedges[i];
816: DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);
818: Yfti = branch->yft[0];
819: Yftr = branch->yft[1];
821: DMNetworkGetConnectedVertices(networkdm,e,&cone);
822: vfrom = cone[0];
823: vto = cone[1];
825: DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom);
826: DMNetworkGetLocalVecOffset(networkdm,vto,0,&offsetto);
line829">829: Vfr = xarr[offse
line830">830: Vfi = xarr[offsetf
line831">831: Vtr = xarr[off
line832">832: Vti = xarr[offse
line834">834: if (vfrom
line835">835: farr[offsetfrom] += Yftr*Vti + Yf
line836">836: farr[offsetfrom+1] += Yftr*Vtr - Yf
line837">837: } else
line838">838: farr[offsetto] += Yftr*Vfi + Yf
line839">839: farr[offsetto+1] += Yftr*Vfr - Yf
line840">840:
line841">841:
line842">842: } else if (key
line843">843: if (!ghost
/* Generator variables */
/* Generator dq axis currents */
line846">846: PetscScalar Vd,Vq,IGr,IGi,Zdq_inv[
/* Generator parameters */
line849">849: gen = (Gen*)(comp
line850">850: DMNetworkGetLocalVecOffset(networkdm,v,j,&offs
/* Generator state variables */
line853">853: Eqp = xarr[offs
line854">854: Edp = xarr[offset
line855">855: delta = xarr[offset
/* w = xarr[idx+3]; not being used */
line857">857: Id = xarr[offset
line858">858: Iq = xarr[offset
/* Generator parameters */
line861">861: Xdp = gen-&
line862">862: Xqp = gen-&
line863">863: Rs = gen-
/* Set generator differential equation residual functions to zero */
line866">866: farr[offsetgen]
line867">867: farr[offsetgen+
line868">868: farr[offsetgen+
line869">869: farr[offsetgen+
line871">871: ri2dq(Vr,Vi,delta,&Vd,&a
/* Algebraic equations for stator currents */
line874">874: det = Rs*Rs + X
line876">876: Zdq_inv[0] =
line877">877: Zdq_inv[1] = X
line878">878: Zdq_inv[2] = -X
line879">879: Zdq_inv[3] =
line881">881: farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq
line882">882: farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq
/* Add generator current injection to network */
line885">885: dq2ri(Id,Iq,delta,&IGr,&am
line887">887: farr[offsetbus]
line888">888: farr[offsetbus+1]
/* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */
line892">892:
line893">893: } else if (key
line894">894: if (!ghost
line895">895: PetscInt off
line896">896: DMNetworkGetLocalVecOffset(networkdm,v,j,&offs
/* Set exciter differential equation residual functions equal to zero*/
line898">898: farr[offsetex
line899">899: farr[offsetexc+
line900">900: farr[offsetexc+
line901">901:
line902">902: } else if (key
line903">903: if (!ghost
line904">904: PetscInt k,ld_nsegsp,ld_
line905">905: PetscScalar *ld_alphap,*ld_betap,*ld_alphaq,*ld_betaq,PD0,QD0,I
line907">907: load = (Load*)(comp
/* Load Parameters */
line910">910: ld_nsegsp = load->ld_
line911">911: ld_alphap = load->ld_
line912">912: ld_betap = load->ld
line913">913: ld_nsegsq = load->ld_
line914">914: ld_alphaq = load->ld_
line915">915: ld_betaq = load->ld
line917">917: PD0 = load-&
line918">918: QD0 = load-&
line920">920: Vm = PetscSqrtScalar(Vr*Vr +
line921">921: Vm2 =
line922">922: Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi
line923">923: PD = QD
line924">924: for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_bet
line925">925: for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_bet
/* Load currents */
line928">928: IDr = (PD*Vr + QD*V
line929">929: IDi = (-QD*Vr + PD*V
line931">931: farr[offsetbus]
line932">932: farr[offsetbus+1]
line933">933:
line934">934:
line935">935:
line936">936:
line938">938: VecRestoreArrayRead(localX,&
line939">939: VecRestoreArray(localF,&
line940">940: DMRestoreLocalVector(networkdm,&l
line942">942: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES<
line943">943: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES<
line944">944: DMRestoreLocalVector(networkdm,&l
line945">945: return
line946">946
line948">948: int main(int argc,char ** argv)
line949">949
line951">951: PetscInt i,j,*edgelist= NULL,eStart,eEnd,vStar
line952">952: PetscInt genj,excj,loadj,component
/* No. of copies (default = 1) */
line954">954: PetscMPIInt siz
line955">955: Vec X
line956">956: TS
line957">957: SNES snes_al
line958">958: Bus
line959">959: Branch *
line960">960: Gen
line961">961: Exc
line962">962: Load
line963">963: DM net
line964">964: #if defined(PETSC_USE_LOG)
line965">965: PetscLogStage
line966">966: #endif
line967">967: Userctx
line968">968: KSP
line969">969: PC
line970">970: PetscInt numEdg
line972">972: PetscInitialize(&argc,&argv,"ex9busnetworkops"
line973">973: PetscOptionsGetInt(NULL,NULL,"-nc",&nc
line974">974: MPI_Comm_size(PETSC_COMM_WORLD,&
line975">975: MPI_Comm_rank(PETSC_COMM_WORLD,&
/* Read initial voltage vector and Ybus */
line978">978: if (rank
line979">979: read_data(nc,&gen,&exc,&load,&bus,&branch,&edg
line980">980:
line982">982: DMNetworkCreate(PETSC_COMM_WORLD,&netw
line983">983: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(Branch),&componentk
line984">984: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(Bus),&componentk
line985">985: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(Gen),&componentk
line986">986: DMNetworkRegisterComponent(networkdm,"excstruct",sizeof(Exc),&componentk
line987">987: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(Load),&componentk
line989">989: PetscLogStageRegister("Create network",&s
line990">990: PetscLogStagePush(s
/* Set local number of edges and edge connectivity */
line993">993: if (rank == 0) numEdges = NBRANCH*nc+
line994">994: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE<
line995">995: DMNetworkAddSubnetwork(networkdm,NULL,numEdges,edgelist
/* Set up the network layout */
line998">998: DMNetworkLayoutSetUp(netw
line1000">1000: if (rank
line1001">1001: PetscFree(edg
line1002">1002:
/* Add network components (physical parameters of nodes and branches) and number of variables */
line1005">1005: if (rank
line1006">1006: DMNetworkGetEdgeRange(networkdm,&eStart,&
line1007">1007: genj=0; loadj=0;
line1008">1008: for (i = eStart; i < eEnd;
line1009">1009: DMNetworkAddComponent(networkdm,i,componentkey[0],&branch[i-eSta
line1010">1010:
line1012">1012: DMNetworkGetVertexRange(networkdm,&vStart,&
line1014">1014: for (i = vStart; i < vEnd;
line1015">1015: DMNetworkAddComponent(networkdm,i,componentkey[1],&bus[i-vSta
line1016">1016: if (bus[i-vStart].no
line1017">1017: for (j = 0; j < bus[i-vStart].nofgen;
/* Add generator */
line1019">1019: DMNetworkAddComponent(networkdm,i,componentkey[2],&gen[genj
/* Add exciter */
line1021">1021: DMNetworkAddComponent(networkdm,i,componentkey[3],&exc[excj
line1022">1022:
line1023">1023:
line1024">1024: if (bus[i-vStart].nof
line1025">1025: for (j=0; j < bus[i-vStart].nofload;
line1026">1026: DMNetworkAddComponent(networkdm,i,componentkey[4],&load[loadj
line1027">1027:
line1028">1028:
line1029">1029:
line1030">1030:
line1032">1032: DMSetUp(netw
line1034">1034: if (rank
line1035">1035: PetscFree5(bus,gen,load,branc
line1036">1036:
/* for parallel options: Network partitioning and distribution of data */
line1039">1039: if (size &g
line1040">1040: DMNetworkDistribute(&networ
line1041">1041:
line1042">1042: PetscLogStagePop
line1044">1044: DMCreateGlobalVector(networkdm,&
line1046">1046: SetInitialGuess(networ
/* Options for fault simulation */
line1049">1049: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options",""
line1050">1050: user.tfaulton
line1051">1051: user.tfaultoff
line1052">1052: user.Rfault =
line1053">1053: user.faultbu
line1054">1054: PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton
line1055">1055: PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff
line1056">1056: PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus
line1057">1057: user.t0
line1058">1058: user.tmax
line1059">1059: PetscOptionsReal("-t0","","",user.t0,&user.t0
line1060">1060: PetscOptionsReal("-tmax","","",user.tmax,&user.tmax
line1062">1062: PetscMalloc1(18*nc,&user.ybus
line1063">1063: for (i = 0; i < 18*nc;
line1064">1064: user.ybusfault[
line1065">1065:
line1066">1066: user.ybusfault[user.faultbus*2+1] = 1/user.
line1067">1067: PetscOptionsEnd
/* Setup TS solver */
/*--------------------------------------------------------*/
line1071">1071: TSCreate(PETSC_COMM_WORLD,&a
line1072">1072: TSSetDM(ts,(DM)netw
line1073">1073: TSSetType(ts,TSC
line1075">1075: TSGetSNES(ts,&
line1076">1076: SNESGetKSP(snes,&am
line1077">1077: KSPGetPC(ksp,&a
line1078">1078: PCSetType(pc,PCBJACOB
line1080">1080: TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1081">1081: TSSetMaxTime(ts,user.tfa
line1082">1082: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1083">1083: TSSetTimeStep(ts
line1084">1084: TSSetFromOptions
/*user.alg_flg = PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix.
eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */
line1088">1088: user.alg_flg = PETSC_FAL
/* Prefault period */
line1091">1091: if (rank
line1092">1092: PetscPrintf(PETSC_COMM_SELF,"... (1) Prefault period ... \n"
line1093">1093:
line1095">1095: TSSetSolution
line1096">1096: TSSetUp
line1097">1097: TSSolve
/* Create the nonlinear solver for solving the algebraic system */
line1100">1100: VecDuplicate(X,&
line1102">1102: SNESCreate(PETSC_COMM_WORLD,&sne
line1103">1103: SNESSetDM(snes_alg,(DM)netw
line1104">1104: SNESSetFunction(snes_alg,F_alg,AlgFunction,&
line1105">1105: SNESSetOptionsPrefix(snes_alg,"alg_"
line1106">1106: SNESSetFromOptions(sne
/* Apply disturbance - resistive fault at user.faultbus */
/* This is done by adding shunt conductance to the diagonal location
in the Ybus matrix */
line1111">1111: user.alg_flg = PETSC_TR
/* Solve the algebraic equations */
line1114">1114: if (rank
line1115">1115: PetscPrintf(PETSC_COMM_SELF,"\n... (2) Apply disturbance, solve algebraic equations ... \n"
line1116">1116:
line1117">1117: SNESSolve(snes_alg,N
/* Disturbance period */
line1120">1120: TSSetTime(ts,user.tfa
line1121">1121: TSSetMaxTime(ts,user.tfau
line1122">1122: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1123">1123: TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1125">1125: user.alg_flg = PETSC_TR
line1126">1126: if (rank
line1127">1127: PetscPrintf(PETSC_COMM_SELF,"\n... (3) Disturbance period ... \n"
line1128">1128:
line1129">1129: TSSolve
/* Remove the fault */
line1132">1132: SNESSetFunction(snes_alg,F_alg,AlgFunction,&
line1134">1134: user.alg_flg = PETSC_FAL
/* Solve the algebraic equations */
line1136">1136: if (rank
line1137">1137: PetscPrintf(PETSC_COMM_SELF,"\n... (4) Remove fault, solve algebraic equations ... \n"
line1138">1138:
line1139">1139: SNESSolve(snes_alg,N
line1140">1140: SNESDestroy(&sne
/* Post-disturbance period */
line1143">1143: TSSetTime(ts,user.tfau
line1144">1144: TSSetMaxTime(ts,user
line1145">1145: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1146">1146: TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1148">1148: user.alg_flg = PETSC_FAL
line1149">1149: if (rank
line1150">1150: PetscPrintf(PETSC_COMM_SELF,"\n... (5) Post-disturbance period ... \n"
line1151">1151:
line1152">1152: TSSolve
line1154">1154: PetscFree(user.ybus
line1155">1155: VecDestroy(&
line1156">1156: VecDestroy(&
line1157">1157: DMDestroy(&netw
line1158">1158: TSDestroy(&a
line1159">1159: PetscFinalize
line1160">1160: return
line1161">1161:
/*TEST
build:
requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
test:
args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
localrunfiles: X.bin Ybus.bin ex9busnetworkops
test:
suffix: 2
nsize: 2
args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
localrunfiles: X.bin Ybus.bin ex9busnetworkops
TEST*/