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: /* T
11: Concepts: DMNetwork
12: Concepts: PETSc TS solver
14: This example was modified from ex9busdmnetwork.c.
15: */
17: #include <petscts.h>
18: #include <petscdmnetwork.h>
20: #define FREQ 60
21: #define W_S (2*PETSC_PI*FREQ)
22: #define NGEN 3 /* No. of generators in the 9 bus system */
23: #define NLOAD 3 /* No. of loads in the 9 bus system */
24: #define NBUS 9 /* No. of buses in the 9 bus system */
25: #define NBRANCH 9 /* No. of branches in the 9 bus system */
27: typedef struct {
28: PetscInt id; /* Bus Number or extended bus name*/
29: PetscScalar mbase; /* MVA base of the machine */
30: PetscScalar PG; /* Generator active power output */
31: PetscScalar QG; /* Generator reactive power output */
33: /* Generator constants */
34: PetscScalar H; /* Inertia constant */
35: PetscScalar Rs; /* Stator Resistance */
36: PetscScalar Xd; /* d-axis reactance */
37: PetscScalar Xdp; /* d-axis transient reactance */
38: PetscScalar Xq; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
39: PetscScalar Xqp; /* q-axis transient reactance */
40: PetscScalar Td0p; /* d-axis open circuit time constant */
41: PetscScalar Tq0p; /* q-axis open circuit time constant */
42: PetscScalar M; /* M = 2*H/W_S */
43: PetscScalar D; /* D = 0.1*M */
44: PetscScalar TM; /* Mechanical Torque */
45: } Gen;
47: typedef struct {
48: /* Exciter system constants */
49: PetscScalar KA ; /* Voltage regulator gain constant */
50: PetscScalar TA; /* Voltage regulator time constant */
51: PetscScalar KE; /* Exciter gain constant */
52: PetscScalar TE; /* Exciter time constant */
53: PetscScalar KF; /* Feedback stabilizer gain constant */
54: PetscScalar TF; /* Feedback stabilizer time constant */
55: PetscScalar k1,k2; /* calculating the saturation function SE = k1*exp(k2*Efd) */
56: PetscScalar Vref; /* Voltage regulator voltage setpoint */
57: } Exc;
59: typedef struct {
60: PetscInt id; /* node id */
61: PetscInt nofgen; /* Number of generators at the bus*/
62: PetscInt nofload; /* Number of load at the bus*/
63: PetscScalar yff[2]; /* yff[0]= imaginary part of admittance, yff[1]=real part of admittance*/
64: PetscScalar vr; /* Real component of bus voltage */
65: PetscScalar vi; /* Imaginary component of bus voltage */
66: } Bus;
68: /* Load constants
69: We use a composite load model that describes the load and reactive powers at each time instant as follows
70: P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
71: Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
72: where
73: id - index of the load
74: ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
75: ld_alphap,ld_alphap - Percentage contribution (weights) or loads
76: P_D0 - Real power load
77: Q_D0 - Reactive power load
78: Vm(t) - Voltage magnitude at time t
79: Vm0 - Voltage magnitude at t = 0
80: ld_betap, ld_betaq - exponents describing the load model for real and reactive part
82: Note: All loads have the same characteristic currently.
83: */
84: typedef struct {
85: PetscInt id; /* bus id */
86: PetscInt ld_nsegsp,ld_nsegsq;
87: PetscScalar PD0, QD0;
88: PetscScalar ld_alphap[3]; /* ld_alphap=[1,0,0], an array, not a value, so use *ld_alphap; */
89: PetscScalar ld_betap[3],ld_alphaq[3],ld_betaq[3];
90: } Load;
92: typedef struct {
93: PetscInt id; /* node id */
94: PetscScalar yft[2]; /* yft[0]= imaginary part of admittance, yft[1]=real part of admittance*/
95: } Branch;
97: typedef struct {
98: PetscReal tfaulton,tfaultoff; /* Fault on and off times */
99: PetscReal t;
100: PetscReal t0,tmax; /* initial time and final time */
101: PetscInt faultbus; /* Fault bus */
102: PetscScalar Rfault; /* Fault resistance (pu) */
103: PetscScalar *ybusfault;
104: PetscBool alg_flg;
105: } Userctx;
107: /* Used to read data into the DMNetwork components */
108: PetscErrorCode read_data(PetscInt nc, Gen **pgen,Exc **pexc, Load **pload,Bus **pbus, Branch **pbranch, PetscInt **pedgelist)
109: {
110: PetscInt i,j,row[1],col[2];
111: PetscInt *edgelist;
112: PetscInt nofgen[9] = {1,1,1,0,0,0,0,0,0}; /* Buses at which generators are incident */
113: PetscInt nofload[9] = {0,0,0,0,1,1,0,1,0}; /* Buses at which loads are incident */
114: const PetscScalar *varr;
115: PetscScalar M[3],D[3];
116: Bus *bus;
117: Branch *branch;
118: Gen *gen;
119: Exc *exc;
120: Load *load;
121: Mat Ybus;
122: Vec V0;
124: /*10 parameters*/
125: /* Generator real and reactive powers (found via loadflow) */
126: static const PetscScalar PG[3] = {0.716786142395021,1.630000000000000,0.850000000000000};
127: static const PetscScalar QG[3] = {0.270702180178785,0.066120127797275,-0.108402221791588};
129: /* Generator constants */
130: static const PetscScalar H[3] = {23.64,6.4,3.01}; /* Inertia constant */
131: static const PetscScalar Rs[3] = {0.0,0.0,0.0}; /* Stator Resistance */
132: static const PetscScalar Xd[3] = {0.146,0.8958,1.3125}; /* d-axis reactance */
133: static const PetscScalar Xdp[3] = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */
134: 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 */
135: static const PetscScalar Xqp[3] = {0.0969,0.1969,0.25}; /* q-axis transient reactance */
136: static const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */
137: static const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */
139: /* Exciter system constants (8 parameters)*/
140: static const PetscScalar KA[3] = {20.0,20.0,20.0}; /* Voltage regulartor gain constant */
141: static const PetscScalar TA[3] = {0.2,0.2,0.2}; /* Voltage regulator time constant */
142: static const PetscScalar KE[3] = {1.0,1.0,1.0}; /* Exciter gain constant */
143: static const PetscScalar TE[3] = {0.314,0.314,0.314}; /* Exciter time constant */
144: static const PetscScalar KF[3] = {0.063,0.063,0.063}; /* Feedback stabilizer gain constant */
145: static const PetscScalar TF[3] = {0.35,0.35,0.35}; /* Feedback stabilizer time constant */
146: static const PetscScalar k1[3] = {0.0039,0.0039,0.0039};
147: static const PetscScalar k2[3] = {1.555,1.555,1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */
149: /* Load constants */
150: static const PetscScalar PD0[3] = {1.25,0.9,1.0};
151: static const PetscScalar QD0[3] = {0.5,0.3,0.35};
152: static const PetscScalar ld_alphaq[3] = {1,0,0};
153: static const PetscScalar ld_betaq[3] = {2,1,0};
154: static const PetscScalar ld_betap[3] = {2,1,0};
155: static const PetscScalar ld_alphap[3] = {1,0,0};
156: PetscInt ld_nsegsp[3] = {3,3,3};
157: PetscInt ld_nsegsq[3] = {3,3,3};
158: PetscViewer Xview,Ybusview;
159: PetscInt neqs_net,m,n;
162: /* Read V0 and Ybus from files */
163: PetscViewerBinaryOpen(PETSC_COMM_SELF,"X.bin",FILE_MODE_READ,&Xview);
164: PetscViewerBinaryOpen(PETSC_COMM_SELF,"Ybus.bin",FILE_MODE_READ,&Ybusview);
165: VecCreate(PETSC_COMM_SELF,&V0);
166: VecLoad(V0,Xview);
168: MatCreate(PETSC_COMM_SELF,&Ybus);
169: MatSetType(Ybus,MATBAIJ);
170: MatLoad(Ybus,Ybusview);
172: /* Destroy unnecessary stuff */
173: PetscViewerDestroy(&Xview);
174: PetscViewerDestroy(&Ybusview);
176: MatGetLocalSize(Ybus,&m,&n);
177: neqs_net = 2*NBUS; /* # eqs. for network subsystem */
180: M[0] = 2*H[0]/W_S;
181: M[1] = 2*H[1]/W_S;
182: M[2] = 2*H[2]/W_S;
183: D[0] = 0.1*M[0];
184: D[1] = 0.1*M[1];
185: D[2] = 0.1*M[2];
187: /* Alocate memory for bus, generators, exciter, loads and branches */
188: PetscCalloc5(NBUS*nc,&bus,NGEN*nc,&gen,NLOAD*nc,&load,NBRANCH*nc+(nc-1),&branch,NGEN*nc,&exc);
190: VecGetArrayRead(V0,&varr);
192: /* read bus data */
193: for (i = 0; i < nc; i++) {
194: for (j = 0; j < NBUS; j++) {
195: bus[i*9+j].id = i*9+j;
196: bus[i*9+j].nofgen = nofgen[j];
197: bus[i*9+j].nofload = nofload[j];
198: bus[i*9+j].vr = varr[2*j];
199: bus[i*9+j].vi = varr[2*j+1];
200: row[0] = 2*j;
201: col[0] = 2*j;
202: col[1] = 2*j+1;
203: /* real and imaginary part of admittance from Ybus into yff */
204: MatGetValues(Ybus,1,row,2,col,bus[i*9+j].yff);
205: }
206: }
208: /* read generator data */
209: for (i = 0; i<nc; i++) {
210: for (j = 0; j < NGEN; j++) {
211: gen[i*3+j].id = i*3+j;
212: gen[i*3+j].PG = PG[j];
213: gen[i*3+j].QG = QG[j];
214: gen[i*3+j].H = H[j];
215: gen[i*3+j].Rs = Rs[j];
216: gen[i*3+j].Xd = Xd[j];
217: gen[i*3+j].Xdp = Xdp[j];
218: gen[i*3+j].Xq = Xq[j];
219: gen[i*3+j].Xqp = Xqp[j];
220: gen[i*3+j].Td0p = Td0p[j];
221: gen[i*3+j].Tq0p = Tq0p[j];
222: gen[i*3+j].M = M[j];
223: gen[i*3+j].D = D[j];
224: }
225: }
227: for (i = 0; i < nc; i++) {
228: for (j = 0; j < NGEN; j++) {
229: /* exciter system */
230: exc[i*3+j].KA = KA[j];
231: exc[i*3+j].TA = TA[j];
232: exc[i*3+j].KE = KE[j];
233: exc[i*3+j].TE = TE[j];
234: exc[i*3+j].KF = KF[j];
235: exc[i*3+j].TF = TF[j];
236: exc[i*3+j].k1 = k1[j];
237: exc[i*3+j].k2 = k2[j];
238: }
239: }
241: /* read load data */
242: for (i = 0; i<nc; i++) {
243: for (j = 0; j < NLOAD; j++) {
244: load[i*3+j].id = i*3+j;
245: load[i*3+j].PD0 = PD0[j];
246: load[i*3+j].QD0 = QD0[j];
247: load[i*3+j].ld_nsegsp = ld_nsegsp[j];
249: load[i*3+j].ld_alphap[0] = ld_alphap[0];
250: load[i*3+j].ld_alphap[1] = ld_alphap[1];
251: load[i*3+j].ld_alphap[2] = ld_alphap[2];
253: load[i*3+j].ld_alphaq[0] = ld_alphaq[0];
254: load[i*3+j].ld_alphaq[1] = ld_alphaq[1];
255: load[i*3+j].ld_alphaq[2] = ld_alphaq[2];
257: load[i*3+j].ld_betap[0] = ld_betap[0];
258: load[i*3+j].ld_betap[1] = ld_betap[1];
259: load[i*3+j].ld_betap[2] = ld_betap[2];
260: load[i*3+j].ld_nsegsq = ld_nsegsq[j];
262: load[i*3+j].ld_betaq[0] = ld_betaq[0];
263: load[i*3+j].ld_betaq[1] = ld_betaq[1];
264: load[i*3+j].ld_betaq[2] = ld_betaq[2];
265: }
266: }
267: PetscCalloc1(2*NBRANCH*nc+2*(nc-1),&edgelist);
269: /* read edgelist */
270: for (i = 0; i<nc; i++) {
271: for (j = 0; j < NBRANCH; j++) {
272: switch (j) {
273: case 0:
274: edgelist[i*18+2*j] = 0+9*i;
275: edgelist[i*18+2*j+1] = 3+9*i;
276: break;
277: case 1:
278: edgelist[i*18+2*j] = 1+9*i;
279: edgelist[i*18+2*j+1] = 6+9*i;
280: break;
281: case 2:
282: edgelist[i*18+2*j] = 2+9*i;
283: edgelist[i*18+2*j+1] = 8+9*i;
284: break;
285: case 3:
286: edgelist[i*18+2*j] = 3+9*i;
287: edgelist[i*18+2*j+1] = 4+9*i;
288: break;
289: case 4:
290: edgelist[i*18+2*j] = 3+9*i;
291: edgelist[i*18+2*j+1] = 5+9*i;
292: break;
293: case 5:
294: edgelist[i*18+2*j] = 4+9*i;
295: edgelist[i*18+2*j+1] = 6+9*i;
296: break;
297: case 6:
298: edgelist[i*18+2*j] = 5+9*i;
299: edgelist[i*18+2*j+1] = 8+9*i;
300: break;
301: case 7:
302: edgelist[i*18+2*j] = 6+9*i;
303: edgelist[i*18+2*j+1] = 7+9*i;
304: break;
305: case 8:
306: edgelist[i*18+2*j] = 7+9*i;
307: edgelist[i*18+2*j+1] = 8+9*i;
308: break;
309: default:
310: break;
311: }
312: }
313: }
315: /* 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 */
316: for (i = 1; i<nc; i++) {
317: edgelist[18*nc+2*(i-1)] = 8+(i-1)*9;
318: edgelist[18*nc+2*(i-1)+1] = 9*i;
320: /* adding admittances to the off-diagonal elements */
321: branch[9*nc+(i-1)].id = 9*nc+(i-1);
322: branch[9*nc+(i-1)].yft[0] = 17.3611;
323: branch[9*nc+(i-1)].yft[1] = -0.0301407;
325: /* subtracting admittances from the diagonal elements */
326: bus[9*i-1].yff[0] -= 17.3611;
327: bus[9*i-1].yff[1] -= -0.0301407;
329: bus[9*i].yff[0] -= 17.3611;
330: bus[9*i].yff[1] -= -0.0301407;
331: }
333: /* read branch data */
334: for (i = 0; i<nc; i++) {
335: for (j = 0; j < NBRANCH; j++) {
336: branch[i*9+j].id = i*9+j;
338: row[0] = edgelist[2*j]*2;
339: col[0] = edgelist[2*j+1]*2;
340: col[1] = edgelist[2*j+1]*2+1;
341: MatGetValues(Ybus,1,row,2,col,branch[i*9+j].yft);/*imaginary part of admittance*/
342: }
343: }
345: *pgen = gen;
346: *pexc = exc;
347: *pload = load;
348: *pbus = bus;
349: *pbranch = branch;
350: *pedgelist = edgelist;
352: VecRestoreArrayRead(V0,&varr);
354: /* Destroy unnecessary stuff */
355: MatDestroy(&Ybus);
356: VecDestroy(&V0);
357: return 0;
358: }
360: PetscErrorCode SetInitialGuess(DM networkdm, Vec X)
361: {
362: Bus *bus;
363: Gen *gen;
364: Exc *exc;
365: PetscInt v,vStart,vEnd,offset;
366: PetscInt key,numComps,j;
367: PetscBool ghostvtex;
368: Vec localX;
369: PetscScalar *xarr;
370: PetscScalar Vr=0,Vi=0,Vm=0,Vm2; /* Terminal voltage variables */
371: PetscScalar IGr, IGi; /* Generator real and imaginary current */
372: PetscScalar Eqp,Edp,delta; /* Generator variables */
373: PetscScalar Efd=0,RF,VR; /* Exciter variables */
374: PetscScalar Vd,Vq; /* Generator dq axis voltages */
375: PetscScalar Id,Iq; /* Generator dq axis currents */
376: PetscScalar theta; /* Generator phase angle */
377: PetscScalar SE;
378: void* component;
380: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
381: DMGetLocalVector(networkdm,&localX);
383: VecSet(X,0.0);
384: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
385: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
387: VecGetArray(localX,&xarr);
389: for (v = vStart; v < vEnd; v++) {
390: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
391: if (ghostvtex) continue;
393: DMNetworkGetNumComponents(networkdm,v,&numComps);
394: for (j=0; j < numComps; j++) {
395: DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL);
396: if (key == 1) {
397: bus = (Bus*)(component);
399: DMNetworkGetLocalVecOffset(networkdm,v,j,&offset);
400: xarr[offset] = bus->vr;
401: xarr[offset+1] = bus->vi;
403: Vr = bus->vr;
404: Vi = bus->vi;
405: } else if (key == 2) {
406: gen = (Gen*)(component);
407: DMNetworkGetLocalVecOffset(networkdm,v,j,&offset);
408: Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi);
409: Vm2 = Vm*Vm;
410: /* Real part of gen current */
411: IGr = (Vr*gen->PG + Vi*gen->QG)/Vm2;
412: /* Imaginary part of gen current */
413: IGi = (Vi*gen->PG - Vr*gen->QG)/Vm2;
415: /* Machine angle */
416: delta = atan2(Vi+gen->Xq*IGr,Vr-gen->Xq*IGi);
417: theta = PETSC_PI/2.0 - delta;
419: /* d-axis stator current */
420: Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta);
422: /* q-axis stator current */
423: Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta);
425: Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta);
426: Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta);
428: /* d-axis transient EMF */
429: Edp = Vd + gen->Rs*Id - gen->Xqp*Iq;
431: /* q-axis transient EMF */
432: Eqp = Vq + gen->Rs*Iq + gen->Xdp*Id;
434: gen->TM = gen->PG;
436: xarr[offset] = Eqp;
437: xarr[offset+1] = Edp;
438: xarr[offset+2] = delta;
439: xarr[offset+3] = W_S;
440: xarr[offset+4] = Id;
441: xarr[offset+5] = Iq;
443: Efd = Eqp + (gen->Xd - gen->Xdp)*Id;
445: } else if (key == 3) {
446: exc = (Exc*)(component);
447: DMNetworkGetLocalVecOffset(networkdm,v,j,&offset);
449: SE = exc->k1*PetscExpScalar(exc->k2*Efd);
450: VR = exc->KE*Efd + SE;
451: RF = exc->KF*Efd/exc->TF;
453: xarr[offset] = Efd;
454: xarr[offset+1] = RF;
455: xarr[offset+2] = VR;
457: exc->Vref = Vm + (VR/exc->KA);
458: }
459: }
460: }
461: VecRestoreArray(localX,&xarr);
462: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
463: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
464: DMRestoreLocalVector(networkdm,&localX);
465: return 0;
466: }
468: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
469: PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr,PetscScalar *Fi)
470: {
471: *Fr = Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta);
472: *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta);
473: return 0;
474: }
476: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
477: PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd,PetscScalar *Fq)
478: {
479: *Fd = Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta);
480: *Fq = Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta);
481: return 0;
482: }
484: /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
485: PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx *user)
486: {
487: DM networkdm;
488: Vec localX,localXdot,localF;
489: PetscInt vfrom,vto,offsetfrom,offsetto;
490: PetscInt v,vStart,vEnd,e;
491: PetscScalar *farr;
492: PetscScalar Vd=0,Vq=0,SE;
493: const PetscScalar *xarr,*xdotarr;
494: void* component;
495: PetscScalar Vr=0, Vi=0;
497: VecSet(F,0.0);
499: TSGetDM(ts,&networkdm);
500: DMGetLocalVector(networkdm,&localF);
501: DMGetLocalVector(networkdm,&localX);
502: DMGetLocalVector(networkdm,&localXdot);
503: VecSet(localF,0.0);
505: /* update ghost values of localX and localXdot */
506: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
507: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
509: DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
510: DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);
512: VecGetArrayRead(localX,&xarr);
513: VecGetArrayRead(localXdot,&xdotarr);
514: VecGetArray(localF,&farr);
516: 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: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
532: DMNetworkGetNumComponents(networkdm,v,&numComps);
534: for (j = 0; j < numComps; j++) {
535: DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL);
536: if (key == 1) {
537: PetscInt nconnedges;
538: const PetscInt *connedges;
540: bus = (Bus*)(component);
541: 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: 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: DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);
574: Yfti = branch->yft[0];
575: Yftr = branch->yft[1];
577: DMNetworkGetConnectedVertices(networkdm,e,&cone);
579: vfrom = cone[0];
580: vto = cone[1];
582: DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom);
583: 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: 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: 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: 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: dq2ri(Id,Iq,delta,&IGr,&IGi);
654: /* Add generator current injection to network */
655: farr[offsetbus] -= IGi;
656: farr[offsetbus+1] -= IGr;
658: }
659: } else if (key == 3) {
660: if (!ghostvtex) {
661: PetscScalar k1,k2,KE,TE,TF,KA,KF,Vref,TA; /* Generator parameters */
662: PetscScalar Efd,RF,VR; /* Exciter variables */
664: exc = (Exc*)(component);
665: DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetexc);
667: Efd = xarr[offsetexc];
668: RF = xarr[offsetexc+1];
669: VR = xarr[offsetexc+2];
671: k1 = exc->k1;
672: k2 = exc->k2;
673: KE = exc->KE;
674: TE = exc->TE;
675: TF = exc->TF;
676: KA = exc->KA;
677: KF = exc->KF;
678: Vref = exc->Vref;
679: TA = exc->TA;
681: Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);
682: SE = k1*PetscExpScalar(k2*Efd);
684: /* Exciter differential equations */
685: farr[offsetexc] = (KE*Efd + SE - VR)/TE + xdotarr[offsetexc];
686: farr[offsetexc+1] = (RF - KF*Efd/TF)/TF + xdotarr[offsetexc+1];
687: farr[offsetexc+2] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[offsetexc+2];
689: }
690: } else if (key ==4) {
691: if (!ghostvtex) {
692: PetscInt k;
693: PetscInt ld_nsegsp;
694: PetscInt ld_nsegsq;
695: PetscScalar *ld_alphap;
696: PetscScalar *ld_betap,*ld_alphaq,*ld_betaq,PD0, QD0, IDr,IDi;
698: load = (Load*)(component);
700: /* Load Parameters */
701: ld_nsegsp = load->ld_nsegsp;
702: ld_alphap = load->ld_alphap;
703: ld_betap = load->ld_betap;
704: ld_nsegsq = load->ld_nsegsq;
705: ld_alphaq = load->ld_alphaq;
706: ld_betaq = load->ld_betaq;
707: PD0 = load->PD0;
708: QD0 = load->QD0;
710: Vr = xarr[offsetbus]; /* Real part of generator terminal voltage */
711: Vi = xarr[offsetbus+1]; /* Imaginary part of the generator terminal voltage */
712: Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi);
713: Vm2 = Vm*Vm;
714: Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0);
715: PD = QD = 0.0;
716: for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]);
717: for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]);
719: /* Load currents */
720: IDr = (PD*Vr + QD*Vi)/Vm2;
721: IDi = (-QD*Vr + PD*Vi)/Vm2;
723: /* Load current contribution to the network */
724: farr[offsetbus] += IDi;
725: farr[offsetbus+1] += IDr;
726: }
727: }
728: }
729: }
731: VecRestoreArrayRead(localX,&xarr);
732: VecRestoreArrayRead(localXdot,&xdotarr);
733: VecRestoreArray(localF,&farr);
734: DMRestoreLocalVector(networkdm,&localX);
735: DMRestoreLocalVector(networkdm,&localXdot);
737: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
738: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
739: DMRestoreLocalVector(networkdm,&localF);
740: return 0;
741: }
743: /* This function is used for solving the algebraic system only during fault on and
744: off times. It computes the entire F and then zeros out the part corresponding to
745: differential equations
746: F = [0;g(y)];
747: */
748: PetscErrorCode AlgFunction (SNES snes, Vec X, Vec F, void *ctx)
749: {
750: DM networkdm;
751: Vec localX,localF;
752: PetscInt vfrom,vto,offsetfrom,offsetto;
753: PetscInt v,vStart,vEnd,e;
754: PetscScalar *farr;
755: Userctx *user=(Userctx*)ctx;
756: const PetscScalar *xarr;
757: void* component;
758: PetscScalar Vr=0,Vi=0;
760: VecSet(F,0.0);
761: SNESGetDM(snes,&networkdm);
762: DMGetLocalVector(networkdm,&localF);
763: DMGetLocalVector(networkdm,&localX);
764: VecSet(localF,0.0);
766: /* update ghost values of locaX and locaXdot */
767: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
768: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
770: VecGetArrayRead(localX,&xarr);
771: VecGetArray(localF,&farr);
773: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
775: for (v=vStart; v < vEnd; v++) {
776: PetscInt i,j,offsetbus,offsetgen,key,numComps;
777: PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0=0, Vi0=0, PD, QD;
778: Bus *bus;
779: Gen *gen;
780: Load *load;
781: PetscBool ghostvtex;
783: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
784: DMNetworkGetNumComponents(networkdm,v,&numComps);
786: for (j = 0; j < numComps; j++) {
787: DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL);
788: if (key == 1) {
789: PetscInt nconnedges;
790: const PetscInt *connedges;
792: bus = (Bus*)(component);
793: DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetbus);
794: if (!ghostvtex) {
795: Vr = xarr[offsetbus];
796: Vi = xarr[offsetbus+1];
798: Yffr = bus->yff[1];
799: Yffi = bus->yff[0];
800: if (user->alg_flg) {
801: Yffr += user->ybusfault[bus->id*2+1];
802: Yffi += user->ybusfault[bus->id*2];
803: }
804: Vr0 = bus->vr;
805: Vi0 = bus->vi;
807: farr[offsetbus] += Yffi*Vr + Yffr*Vi;
808: farr[offsetbus+1] += Yffr*Vr - Yffi*Vi;
809: }
810: DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
812: for (i=0; i < nconnedges; i++) {
813: Branch *branch;
814: PetscInt keye;
815: PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
816: const PetscInt *cone;
818: e = connedges[i];
819: DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);
821: Yfti = branch->yft[0];
822: Yftr = branch->yft[1];
824: DMNetworkGetConnectedVertices(networkdm,e,&cone);
825: vfrom = cone[0];
826: vto = cone[1];
828: DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom);
829: DMNetworkGetLocalVecOffset(networkdm,vto,0,&offsetto);
line832">832: Vfr = xarr[offse
line833">833: Vfi = xarr[offsetf
line834">834: Vtr = xarr[off
line835">835: Vti = xarr[offse
line837">837: if (vfrom
line838">838: farr[offsetfrom] += Yftr*Vti + Yf
line839">839: farr[offsetfrom+1] += Yftr*Vtr - Yf
line840">840: } else
line841">841: farr[offsetto] += Yftr*Vfi + Yf
line842">842: farr[offsetto+1] += Yftr*Vfr - Yf
line843">843:
line844">844:
line845">845: } else if (key
line846">846: if (!ghost
/* Generator variables */
/* Generator dq axis currents */
line849">849: PetscScalar Vd,Vq,IGr,IGi,Zdq_inv[
/* Generator parameters */
line852">852: gen = (Gen*)(comp
line853">853: DMNetworkGetLocalVecOffset(networkdm,v,j,&offs
/* Generator state variables */
line856">856: Eqp = xarr[offs
line857">857: Edp = xarr[offset
line858">858: delta = xarr[offset
/* w = xarr[idx+3]; not being used */
line860">860: Id = xarr[offset
line861">861: Iq = xarr[offset
/* Generator parameters */
line864">864: Xdp = gen-&
line865">865: Xqp = gen-&
line866">866: Rs = gen-
/* Set generator differential equation residual functions to zero */
line869">869: farr[offsetgen]
line870">870: farr[offsetgen+
line871">871: farr[offsetgen+
line872">872: farr[offsetgen+
line874">874: ri2dq(Vr,Vi,delta,&Vd,&a
/* Algebraic equations for stator currents */
line877">877: det = Rs*Rs + X
line879">879: Zdq_inv[0] =
line880">880: Zdq_inv[1] = X
line881">881: Zdq_inv[2] = -X
line882">882: Zdq_inv[3] =
line884">884: farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq
line885">885: farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq
/* Add generator current injection to network */
line888">888: dq2ri(Id,Iq,delta,&IGr,&am
line890">890: farr[offsetbus]
line891">891: farr[offsetbus+1]
/* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */
line895">895:
line896">896: } else if (key
line897">897: if (!ghost
line898">898: PetscInt off
line899">899: DMNetworkGetLocalVecOffset(networkdm,v,j,&offs
/* Set exciter differential equation residual functions equal to zero*/
line901">901: farr[offsetex
line902">902: farr[offsetexc+
line903">903: farr[offsetexc+
line904">904:
line905">905: } else if (key
line906">906: if (!ghost
line907">907: PetscInt k,ld_nsegsp,ld_
line908">908: PetscScalar *ld_alphap,*ld_betap,*ld_alphaq,*ld_betaq,PD0,QD0,I
line910">910: load = (Load*)(comp
/* Load Parameters */
line913">913: ld_nsegsp = load->ld_
line914">914: ld_alphap = load->ld_
line915">915: ld_betap = load->ld
line916">916: ld_nsegsq = load->ld_
line917">917: ld_alphaq = load->ld_
line918">918: ld_betaq = load->ld
line920">920: PD0 = load-&
line921">921: QD0 = load-&
line923">923: Vm = PetscSqrtScalar(Vr*Vr +
line924">924: Vm2 =
line925">925: Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi
line926">926: PD = QD
line927">927: for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_bet
line928">928: for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_bet
/* Load currents */
line931">931: IDr = (PD*Vr + QD*V
line932">932: IDi = (-QD*Vr + PD*V
line934">934: farr[offsetbus]
line935">935: farr[offsetbus+1]
line936">936:
line937">937:
line938">938:
line939">939:
line941">941: VecRestoreArrayRead(localX,&
line942">942: VecRestoreArray(localF,&
line943">943: DMRestoreLocalVector(networkdm,&l
line945">945: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES<
line946">946: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES<
line947">947: DMRestoreLocalVector(networkdm,&l
line948">948: return
line949">949
line951">951: int main(int argc,char ** argv)
line952">952
line954">954: PetscInt i,j,*edgelist= NULL,eStart,eEnd,vStar
line955">955: PetscInt genj,excj,loadj,component
/* No. of copies (default = 1) */
line957">957: PetscMPIInt siz
line958">958: Vec X
line959">959: TS
line960">960: SNES snes_al
line961">961: Bus
line962">962: Branch *
line963">963: Gen
line964">964: Exc
line965">965: Load
line966">966: DM net
line967">967: #if defined(PETSC_USE_LOG)
line968">968: PetscLogStage
line969">969: #endif
line970">970: Userctx
line971">971: KSP
line972">972: PC
line973">973: PetscInt numEdg
line975">975: PetscInitialize(&argc,&argv,"ex9busnetworkops"
line976">976: PetscOptionsGetInt(NULL,NULL,"-nc",&nc
line977">977: MPI_Comm_size(PETSC_COMM_WORLD,&
line978">978: MPI_Comm_rank(PETSC_COMM_WORLD,&
/* Read initial voltage vector and Ybus */
line981">981: if (rank
line982">982: read_data(nc,&gen,&exc,&load,&bus,&branch,&edg
line983">983:
line985">985: DMNetworkCreate(PETSC_COMM_WORLD,&netw
line986">986: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(Branch),&componentk
line987">987: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(Bus),&componentk
line988">988: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(Gen),&componentk
line989">989: DMNetworkRegisterComponent(networkdm,"excstruct",sizeof(Exc),&componentk
line990">990: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(Load),&componentk
line992">992: PetscLogStageRegister("Create network",&s
line993">993: PetscLogStagePush(s
/* Set local number of edges and edge connectivity */
line996">996: if (rank == 0) numEdges = NBRANCH*nc+
line997">997: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE<
line998">998: DMNetworkAddSubnetwork(networkdm,NULL,numEdges,edgelist
/* Set up the network layout */
line1001">1001: DMNetworkLayoutSetUp(netw
line1003">1003: if (rank
line1004">1004: PetscFree(edg
line1005">1005:
/* Add network components (physical parameters of nodes and branches) and number of variables */
line1008">1008: if (rank
line1009">1009: DMNetworkGetEdgeRange(networkdm,&eStart,&
line1010">1010: genj=0; loadj=0;
line1011">1011: for (i = eStart; i < eEnd;
line1012">1012: DMNetworkAddComponent(networkdm,i,componentkey[0],&branch[i-eSta
line1013">1013:
line1015">1015: DMNetworkGetVertexRange(networkdm,&vStart,&
line1017">1017: for (i = vStart; i < vEnd;
line1018">1018: DMNetworkAddComponent(networkdm,i,componentkey[1],&bus[i-vSta
line1019">1019: if (bus[i-vStart].no
line1020">1020: for (j = 0; j < bus[i-vStart].nofgen;
/* Add generator */
line1022">1022: DMNetworkAddComponent(networkdm,i,componentkey[2],&gen[genj
/* Add exciter */
line1024">1024: DMNetworkAddComponent(networkdm,i,componentkey[3],&exc[excj
line1025">1025:
line1026">1026:
line1027">1027: if (bus[i-vStart].nof
line1028">1028: for (j=0; j < bus[i-vStart].nofload;
line1029">1029: DMNetworkAddComponent(networkdm,i,componentkey[4],&load[loadj
line1030">1030:
line1031">1031:
line1032">1032:
line1033">1033:
line1035">1035: DMSetUp(netw
line1037">1037: if (rank
line1038">1038: PetscFree5(bus,gen,load,branc
line1039">1039:
/* for parallel options: Network partitioning and distribution of data */
line1042">1042: if (size &g
line1043">1043: DMNetworkDistribute(&networ
line1044">1044:
line1045">1045: PetscLogStagePop
line1047">1047: DMCreateGlobalVector(networkdm,&
line1049">1049: SetInitialGuess(networ
/* Options for fault simulation */
line1052">1052: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options",""
line1053">1053: user.tfaulton
line1054">1054: user.tfaultoff
line1055">1055: user.Rfault =
line1056">1056: user.faultbu
line1057">1057: PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton
line1058">1058: PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff
line1059">1059: PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus
line1060">1060: user.t0
line1061">1061: user.tmax
line1062">1062: PetscOptionsReal("-t0","","",user.t0,&user.t0
line1063">1063: PetscOptionsReal("-tmax","","",user.tmax,&user.tmax
line1065">1065: PetscMalloc1(18*nc,&user.ybus
line1066">1066: for (i = 0; i < 18*nc;
line1067">1067: user.ybusfault[
line1068">1068:
line1069">1069: user.ybusfault[user.faultbus*2+1] = 1/user.
line1070">1070: PetscOptionsEnd
/* Setup TS solver */
/*--------------------------------------------------------*/
line1074">1074: TSCreate(PETSC_COMM_WORLD,&a
line1075">1075: TSSetDM(ts,(DM)netw
line1076">1076: TSSetType(ts,TSC
line1078">1078: TSGetSNES(ts,&
line1079">1079: SNESGetKSP(snes,&am
line1080">1080: KSPGetPC(ksp,&a
line1081">1081: PCSetType(pc,PCBJACOB
line1083">1083: TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1084">1084: TSSetMaxTime(ts,user.tfa
line1085">1085: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1086">1086: TSSetTimeStep(ts
line1087">1087: 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. */
line1091">1091: user.alg_flg = PETSC_FAL
/* Prefault period */
line1094">1094: if (rank
line1095">1095: PetscPrintf(PETSC_COMM_SELF,"... (1) Prefault period ... \n"
line1096">1096:
line1098">1098: TSSetSolution
line1099">1099: TSSetUp
line1100">1100: TSSolve
/* Create the nonlinear solver for solving the algebraic system */
line1103">1103: VecDuplicate(X,&
line1105">1105: SNESCreate(PETSC_COMM_WORLD,&sne
line1106">1106: SNESSetDM(snes_alg,(DM)netw
line1107">1107: SNESSetFunction(snes_alg,F_alg,AlgFunction,&
line1108">1108: SNESSetOptionsPrefix(snes_alg,"alg_"
line1109">1109: SNESSetFromOptions(sne
/* Apply disturbance - resistive fault at user.faultbus */
/* This is done by adding shunt conductance to the diagonal location
in the Ybus matrix */
line1114">1114: user.alg_flg = PETSC_TR
/* Solve the algebraic equations */
line1117">1117: if (rank
line1118">1118: PetscPrintf(PETSC_COMM_SELF,"\n... (2) Apply disturbance, solve algebraic equations ... \n"
line1119">1119:
line1120">1120: SNESSolve(snes_alg,N
/* Disturbance period */
line1123">1123: TSSetTime(ts,user.tfa
line1124">1124: TSSetMaxTime(ts,user.tfau
line1125">1125: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1126">1126: TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1128">1128: user.alg_flg = PETSC_TR
line1129">1129: if (rank
line1130">1130: PetscPrintf(PETSC_COMM_SELF,"\n... (3) Disturbance period ... \n"
line1131">1131:
line1132">1132: TSSolve
/* Remove the fault */
line1135">1135: SNESSetFunction(snes_alg,F_alg,AlgFunction,&
line1137">1137: user.alg_flg = PETSC_FAL
/* Solve the algebraic equations */
line1139">1139: if (rank
line1140">1140: PetscPrintf(PETSC_COMM_SELF,"\n... (4) Remove fault, solve algebraic equations ... \n"
line1141">1141:
line1142">1142: SNESSolve(snes_alg,N
line1143">1143: SNESDestroy(&sne
/* Post-disturbance period */
line1146">1146: TSSetTime(ts,user.tfau
line1147">1147: TSSetMaxTime(ts,user
line1148">1148: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1149">1149: TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1151">1151: user.alg_flg = PETSC_FAL
line1152">1152: if (rank
line1153">1153: PetscPrintf(PETSC_COMM_SELF,"\n... (5) Post-disturbance period ... \n"
line1154">1154:
line1155">1155: TSSolve
line1157">1157: PetscFree(user.ybus
line1158">1158: VecDestroy(&
line1159">1159: VecDestroy(&
line1160">1160: DMDestroy(&netw
line1161">1161: TSDestroy(&a
line1162">1162: PetscFinalize
line1163">1163: return
line1164">1164:
/*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*/