Actual source code: ex9bus.c
2: static char help[] = "Power grid stability analysis of WECC 9 bus system.\n\
3: This example is based on the 9-bus (node) example given in the book Power\n\
4: Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\
5: The power grid in this example consists of 9 buses (nodes), 3 generators,\n\
6: 3 loads, and 9 transmission lines. The network equations are written\n\
7: in current balance form using rectangular coordinates.\n\n";
9: /*
10: The equations for the stability analysis are described by the DAE
12: \dot{x} = f(x,y,t)
13: 0 = g(x,y,t)
15: where the generators are described by differential equations, while the algebraic
16: constraints define the network equations.
18: The generators are modeled with a 4th order differential equation describing the electrical
19: and mechanical dynamics. Each generator also has an exciter system modeled by 3rd order
20: diff. eqns. describing the exciter, voltage regulator, and the feedback stabilizer
21: mechanism.
23: The network equations are described by nodal current balance equations.
24: I(x,y) - Y*V = 0
26: where:
27: I(x,y) is the current injected from generators and loads.
28: Y is the admittance matrix, and
29: V is the voltage vector
30: */
32: /*
33: Include "petscts.h" so that we can use TS solvers. Note that this
34: file automatically includes:
35: petscsys.h - base PETSc routines petscvec.h - vectors
36: petscmat.h - matrices
37: petscis.h - index sets petscksp.h - Krylov subspace methods
38: petscviewer.h - viewers petscpc.h - preconditioners
39: petscksp.h - linear solvers
40: */
42: #include <petscts.h>
43: #include <petscdm.h>
44: #include <petscdmda.h>
45: #include <petscdmcomposite.h>
47: #define freq 60
48: #define w_s (2 * PETSC_PI * freq)
50: /* Sizes and indices */
51: const PetscInt nbus = 9; /* Number of network buses */
52: const PetscInt ngen = 3; /* Number of generators */
53: const PetscInt nload = 3; /* Number of loads */
54: const PetscInt gbus[3] = {0, 1, 2}; /* Buses at which generators are incident */
55: const PetscInt lbus[3] = {4, 5, 7}; /* Buses at which loads are incident */
57: /* Generator real and reactive powers (found via loadflow) */
58: const PetscScalar PG[3] = {0.716786142395021, 1.630000000000000, 0.850000000000000};
59: const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};
60: /* Generator constants */
61: const PetscScalar H[3] = {23.64, 6.4, 3.01}; /* Inertia constant */
62: const PetscScalar Rs[3] = {0.0, 0.0, 0.0}; /* Stator Resistance */
63: const PetscScalar Xd[3] = {0.146, 0.8958, 1.3125}; /* d-axis reactance */
64: const PetscScalar Xdp[3] = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
65: 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 */
66: const PetscScalar Xqp[3] = {0.0969, 0.1969, 0.25}; /* q-axis transient reactance */
67: const PetscScalar Td0p[3] = {8.96, 6.0, 5.89}; /* d-axis open circuit time constant */
68: const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6}; /* q-axis open circuit time constant */
69: PetscScalar M[3]; /* M = 2*H/w_s */
70: PetscScalar D[3]; /* D = 0.1*M */
72: PetscScalar TM[3]; /* Mechanical Torque */
73: /* Exciter system constants */
74: const PetscScalar KA[3] = {20.0, 20.0, 20.0}; /* Voltage regulartor gain constant */
75: const PetscScalar TA[3] = {0.2, 0.2, 0.2}; /* Voltage regulator time constant */
76: const PetscScalar KE[3] = {1.0, 1.0, 1.0}; /* Exciter gain constant */
77: const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
78: const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
79: const PetscScalar TF[3] = {0.35, 0.35, 0.35}; /* Feedback stabilizer time constant */
80: const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
81: const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */
82: const PetscScalar VRMIN[3] = {-4.0, -4.0, -4.0};
83: const PetscScalar VRMAX[3] = {7.0, 7.0, 7.0};
84: PetscInt VRatmin[3];
85: PetscInt VRatmax[3];
87: PetscScalar Vref[3];
88: /* Load constants
89: We use a composite load model that describes the load and reactive powers at each time instant as follows
90: P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
91: Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
92: where
93: ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
94: ld_alphap,ld_alphap - Percentage contribution (weights) or loads
95: P_D0 - Real power load
96: Q_D0 - Reactive power load
97: V_m(t) - Voltage magnitude at time t
98: V_m0 - Voltage magnitude at t = 0
99: ld_betap, ld_betaq - exponents describing the load model for real and reactive part
101: Note: All loads have the same characteristic currently.
102: */
103: const PetscScalar PD0[3] = {1.25, 0.9, 1.0};
104: const PetscScalar QD0[3] = {0.5, 0.3, 0.35};
105: const PetscInt ld_nsegsp[3] = {3, 3, 3};
106: const PetscScalar ld_alphap[3] = {1.0, 0.0, 0.0};
107: const PetscScalar ld_betap[3] = {2.0, 1.0, 0.0};
108: const PetscInt ld_nsegsq[3] = {3, 3, 3};
109: const PetscScalar ld_alphaq[3] = {1.0, 0.0, 0.0};
110: const PetscScalar ld_betaq[3] = {2.0, 1.0, 0.0};
112: typedef struct {
113: DM dmgen, dmnet; /* DMs to manage generator and network subsystem */
114: DM dmpgrid; /* Composite DM to manage the entire power grid */
115: Mat Ybus; /* Network admittance matrix */
116: Vec V0; /* Initial voltage vector (Power flow solution) */
117: PetscReal tfaulton, tfaultoff; /* Fault on and off times */
118: PetscInt faultbus; /* Fault bus */
119: PetscScalar Rfault;
120: PetscReal t0, tmax;
121: PetscInt neqs_gen, neqs_net, neqs_pgrid;
122: Mat Sol; /* Matrix to save solution at each time step */
123: PetscInt stepnum;
124: PetscReal t;
125: SNES snes_alg;
126: IS is_diff; /* indices for differential equations */
127: IS is_alg; /* indices for algebraic equations */
128: PetscBool setisdiff; /* TS computes truncation error based only on the differential variables */
129: PetscBool semiexplicit; /* If the flag is set then a semi-explicit method is used using TSRK */
130: } Userctx;
132: /*
133: The first two events are for fault on and off, respectively. The following events are
134: to check the min/max limits on the state variable VR. A non windup limiter is used for
135: the VR limits.
136: */
137: PetscErrorCode EventFunction(TS ts, PetscReal t, Vec X, PetscScalar *fvalue, void *ctx)
138: {
139: Userctx *user = (Userctx *)ctx;
140: Vec Xgen, Xnet;
141: PetscInt i, idx = 0;
142: const PetscScalar *xgen, *xnet;
143: PetscScalar Efd, RF, VR, Vr, Vi, Vm;
146: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
147: DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);
149: VecGetArrayRead(Xgen, &xgen);
150: VecGetArrayRead(Xnet, &xnet);
152: /* Event for fault-on time */
153: fvalue[0] = t - user->tfaulton;
154: /* Event for fault-off time */
155: fvalue[1] = t - user->tfaultoff;
157: for (i = 0; i < ngen; i++) {
158: Efd = xgen[idx + 6];
159: RF = xgen[idx + 7];
160: VR = xgen[idx + 8];
162: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
163: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
164: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
166: if (!VRatmax[i]) {
167: fvalue[2 + 2 * i] = VRMAX[i] - VR;
168: } else {
169: fvalue[2 + 2 * i] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
170: }
171: if (!VRatmin[i]) {
172: fvalue[2 + 2 * i + 1] = VRMIN[i] - VR;
173: } else {
174: fvalue[2 + 2 * i + 1] = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
175: }
176: idx = idx + 9;
177: }
178: VecRestoreArrayRead(Xgen, &xgen);
179: VecRestoreArrayRead(Xnet, &xnet);
181: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
183: return 0;
184: }
186: PetscErrorCode PostEventFunction(TS ts, PetscInt nevents, PetscInt event_list[], PetscReal t, Vec X, PetscBool forwardsolve, void *ctx)
187: {
188: Userctx *user = (Userctx *)ctx;
189: Vec Xgen, Xnet;
190: PetscScalar *xgen, *xnet;
191: PetscInt row_loc, col_loc;
192: PetscScalar val;
193: PetscInt i, idx = 0, event_num;
194: PetscScalar fvalue;
195: PetscScalar Efd, RF, VR;
196: PetscScalar Vr, Vi, Vm;
199: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
200: DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);
202: VecGetArray(Xgen, &xgen);
203: VecGetArray(Xnet, &xnet);
205: for (i = 0; i < nevents; i++) {
206: if (event_list[i] == 0) {
207: /* Apply disturbance - resistive fault at user->faultbus */
208: /* This is done by adding shunt conductance to the diagonal location
209: in the Ybus matrix */
210: row_loc = 2 * user->faultbus;
211: col_loc = 2 * user->faultbus + 1; /* Location for G */
212: val = 1 / user->Rfault;
213: MatSetValues(user->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
214: row_loc = 2 * user->faultbus + 1;
215: col_loc = 2 * user->faultbus; /* Location for G */
216: val = 1 / user->Rfault;
217: MatSetValues(user->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
219: MatAssemblyBegin(user->Ybus, MAT_FINAL_ASSEMBLY);
220: MatAssemblyEnd(user->Ybus, MAT_FINAL_ASSEMBLY);
222: /* Solve the algebraic equations */
223: SNESSolve(user->snes_alg, NULL, X);
224: } else if (event_list[i] == 1) {
225: /* Remove the fault */
226: row_loc = 2 * user->faultbus;
227: col_loc = 2 * user->faultbus + 1;
228: val = -1 / user->Rfault;
229: MatSetValues(user->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
230: row_loc = 2 * user->faultbus + 1;
231: col_loc = 2 * user->faultbus;
232: val = -1 / user->Rfault;
233: MatSetValues(user->Ybus, 1, &row_loc, 1, &col_loc, &val, ADD_VALUES);
235: MatAssemblyBegin(user->Ybus, MAT_FINAL_ASSEMBLY);
236: MatAssemblyEnd(user->Ybus, MAT_FINAL_ASSEMBLY);
238: /* Solve the algebraic equations */
239: SNESSolve(user->snes_alg, NULL, X);
241: /* Check the VR derivatives and reset flags if needed */
242: for (i = 0; i < ngen; i++) {
243: Efd = xgen[idx + 6];
244: RF = xgen[idx + 7];
245: VR = xgen[idx + 8];
247: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
248: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
249: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
251: if (VRatmax[i]) {
252: fvalue = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
253: if (fvalue < 0) {
254: VRatmax[i] = 0;
255: PetscPrintf(PETSC_COMM_SELF, "VR[%" PetscInt_FMT "]: dVR_dt went negative on fault clearing at time %g\n", i, (double)t);
256: }
257: }
258: if (VRatmin[i]) {
259: fvalue = (VR - KA[i] * RF + KA[i] * KF[i] * Efd / TF[i] - KA[i] * (Vref[i] - Vm)) / TA[i];
261: if (fvalue > 0) {
262: VRatmin[i] = 0;
263: PetscPrintf(PETSC_COMM_SELF, "VR[%" PetscInt_FMT "]: dVR_dt went positive on fault clearing at time %g\n", i, (double)t);
264: }
265: }
266: idx = idx + 9;
267: }
268: } else {
269: idx = (event_list[i] - 2) / 2;
270: event_num = (event_list[i] - 2) % 2;
271: if (event_num == 0) { /* Max VR */
272: if (!VRatmax[idx]) {
273: VRatmax[idx] = 1;
274: PetscPrintf(PETSC_COMM_SELF, "VR[%" PetscInt_FMT "]: hit upper limit at time %g\n", idx, (double)t);
275: } else {
276: VRatmax[idx] = 0;
277: PetscPrintf(PETSC_COMM_SELF, "VR[%" PetscInt_FMT "]: freeing variable as dVR_dt is negative at time %g\n", idx, (double)t);
278: }
279: } else {
280: if (!VRatmin[idx]) {
281: VRatmin[idx] = 1;
282: PetscPrintf(PETSC_COMM_SELF, "VR[%" PetscInt_FMT "]: hit lower limit at time %g\n", idx, (double)t);
283: } else {
284: VRatmin[idx] = 0;
285: PetscPrintf(PETSC_COMM_SELF, "VR[%" PetscInt_FMT "]: freeing variable as dVR_dt is positive at time %g\n", idx, (double)t);
286: }
287: }
288: }
289: }
291: VecRestoreArray(Xgen, &xgen);
292: VecRestoreArray(Xnet, &xnet);
294: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
296: return 0;
297: }
299: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
300: PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
301: {
302: *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
303: *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
304: return 0;
305: }
307: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
308: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
309: {
310: *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
311: *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
312: return 0;
313: }
315: /* Saves the solution at each time to a matrix */
316: PetscErrorCode SaveSolution(TS ts)
317: {
318: Userctx *user;
319: Vec X;
320: const PetscScalar *x;
321: PetscScalar *mat;
322: PetscInt idx;
323: PetscReal t;
325: TSGetApplicationContext(ts, &user);
326: TSGetTime(ts, &t);
327: TSGetSolution(ts, &X);
328: idx = user->stepnum * (user->neqs_pgrid + 1);
329: MatDenseGetArray(user->Sol, &mat);
330: VecGetArrayRead(X, &x);
331: mat[idx] = t;
332: PetscArraycpy(mat + idx + 1, x, user->neqs_pgrid);
333: MatDenseRestoreArray(user->Sol, &mat);
334: VecRestoreArrayRead(X, &x);
335: user->stepnum++;
336: return 0;
337: }
339: PetscErrorCode SetInitialGuess(Vec X, Userctx *user)
340: {
341: Vec Xgen, Xnet;
342: PetscScalar *xgen;
343: const PetscScalar *xnet;
344: PetscInt i, idx = 0;
345: PetscScalar Vr, Vi, IGr, IGi, Vm, Vm2;
346: PetscScalar Eqp, Edp, delta;
347: PetscScalar Efd, RF, VR; /* Exciter variables */
348: PetscScalar Id, Iq; /* Generator dq axis currents */
349: PetscScalar theta, Vd, Vq, SE;
351: M[0] = 2 * H[0] / w_s;
352: M[1] = 2 * H[1] / w_s;
353: M[2] = 2 * H[2] / w_s;
354: D[0] = 0.1 * M[0];
355: D[1] = 0.1 * M[1];
356: D[2] = 0.1 * M[2];
358: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
360: /* Network subsystem initialization */
361: VecCopy(user->V0, Xnet);
363: /* Generator subsystem initialization */
364: VecGetArrayWrite(Xgen, &xgen);
365: VecGetArrayRead(Xnet, &xnet);
367: for (i = 0; i < ngen; i++) {
368: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
369: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
370: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
371: Vm2 = Vm * Vm;
372: IGr = (Vr * PG[i] + Vi * QG[i]) / Vm2;
373: IGi = (Vi * PG[i] - Vr * QG[i]) / Vm2;
375: delta = PetscAtan2Real(Vi + Xq[i] * IGr, Vr - Xq[i] * IGi); /* Machine angle */
377: theta = PETSC_PI / 2.0 - delta;
379: Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta); /* d-axis stator current */
380: Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta); /* q-axis stator current */
382: Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
383: Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);
385: Edp = Vd + Rs[i] * Id - Xqp[i] * Iq; /* d-axis transient EMF */
386: Eqp = Vq + Rs[i] * Iq + Xdp[i] * Id; /* q-axis transient EMF */
388: TM[i] = PG[i];
390: /* The generator variables are ordered as [Eqp,Edp,delta,w,Id,Iq] */
391: xgen[idx] = Eqp;
392: xgen[idx + 1] = Edp;
393: xgen[idx + 2] = delta;
394: xgen[idx + 3] = w_s;
396: idx = idx + 4;
398: xgen[idx] = Id;
399: xgen[idx + 1] = Iq;
401: idx = idx + 2;
403: /* Exciter */
404: Efd = Eqp + (Xd[i] - Xdp[i]) * Id;
405: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
406: VR = KE[i] * Efd + SE;
407: RF = KF[i] * Efd / TF[i];
409: xgen[idx] = Efd;
410: xgen[idx + 1] = RF;
411: xgen[idx + 2] = VR;
413: Vref[i] = Vm + (VR / KA[i]);
415: VRatmin[i] = VRatmax[i] = 0;
417: idx = idx + 3;
418: }
419: VecRestoreArrayWrite(Xgen, &xgen);
420: VecRestoreArrayRead(Xnet, &xnet);
422: /* VecView(Xgen,0); */
423: DMCompositeGather(user->dmpgrid, INSERT_VALUES, X, Xgen, Xnet);
424: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
425: return 0;
426: }
428: /* Computes F = [f(x,y);g(x,y)] */
429: PetscErrorCode ResidualFunction(Vec X, Vec F, Userctx *user)
430: {
431: Vec Xgen, Xnet, Fgen, Fnet;
432: const PetscScalar *xgen, *xnet;
433: PetscScalar *fgen, *fnet;
434: PetscInt i, idx = 0;
435: PetscScalar Vr, Vi, Vm, Vm2;
436: PetscScalar Eqp, Edp, delta, w; /* Generator variables */
437: PetscScalar Efd, RF, VR; /* Exciter variables */
438: PetscScalar Id, Iq; /* Generator dq axis currents */
439: PetscScalar Vd, Vq, SE;
440: PetscScalar IGr, IGi, IDr, IDi;
441: PetscScalar Zdq_inv[4], det;
442: PetscScalar PD, QD, Vm0, *v0;
443: PetscInt k;
445: VecZeroEntries(F);
446: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
447: DMCompositeGetLocalVectors(user->dmpgrid, &Fgen, &Fnet);
448: DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);
449: DMCompositeScatter(user->dmpgrid, F, Fgen, Fnet);
451: /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
452: The generator current injection, IG, and load current injection, ID are added later
453: */
454: /* Note that the values in Ybus are stored assuming the imaginary current balance
455: equation is ordered first followed by real current balance equation for each bus.
456: Thus imaginary current contribution goes in location 2*i, and
457: real current contribution in 2*i+1
458: */
459: MatMult(user->Ybus, Xnet, Fnet);
461: VecGetArrayRead(Xgen, &xgen);
462: VecGetArrayRead(Xnet, &xnet);
463: VecGetArrayWrite(Fgen, &fgen);
464: VecGetArrayWrite(Fnet, &fnet);
466: /* Generator subsystem */
467: for (i = 0; i < ngen; i++) {
468: Eqp = xgen[idx];
469: Edp = xgen[idx + 1];
470: delta = xgen[idx + 2];
471: w = xgen[idx + 3];
472: Id = xgen[idx + 4];
473: Iq = xgen[idx + 5];
474: Efd = xgen[idx + 6];
475: RF = xgen[idx + 7];
476: VR = xgen[idx + 8];
478: /* Generator differential equations */
479: fgen[idx] = (-Eqp - (Xd[i] - Xdp[i]) * Id + Efd) / Td0p[i];
480: fgen[idx + 1] = (-Edp + (Xq[i] - Xqp[i]) * Iq) / Tq0p[i];
481: fgen[idx + 2] = w - w_s;
482: fgen[idx + 3] = (TM[i] - Edp * Id - Eqp * Iq - (Xqp[i] - Xdp[i]) * Id * Iq - D[i] * (w - w_s)) / M[i];
484: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
485: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
487: ri2dq(Vr, Vi, delta, &Vd, &Vq);
488: /* Algebraic equations for stator currents */
489: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
491: Zdq_inv[0] = Rs[i] / det;
492: Zdq_inv[1] = Xqp[i] / det;
493: Zdq_inv[2] = -Xdp[i] / det;
494: Zdq_inv[3] = Rs[i] / det;
496: fgen[idx + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
497: fgen[idx + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;
499: /* Add generator current injection to network */
500: dq2ri(Id, Iq, delta, &IGr, &IGi);
502: fnet[2 * gbus[i]] -= IGi;
503: fnet[2 * gbus[i] + 1] -= IGr;
505: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
507: SE = k1[i] * PetscExpScalar(k2[i] * Efd);
509: /* Exciter differential equations */
510: fgen[idx + 6] = (-KE[i] * Efd - SE + VR) / TE[i];
511: fgen[idx + 7] = (-RF + KF[i] * Efd / TF[i]) / TF[i];
512: if (VRatmax[i]) fgen[idx + 8] = VR - VRMAX[i];
513: else if (VRatmin[i]) fgen[idx + 8] = VRMIN[i] - VR;
514: else fgen[idx + 8] = (-VR + KA[i] * RF - KA[i] * KF[i] * Efd / TF[i] + KA[i] * (Vref[i] - Vm)) / TA[i];
516: idx = idx + 9;
517: }
519: VecGetArray(user->V0, &v0);
520: for (i = 0; i < nload; i++) {
521: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
522: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
523: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
524: Vm2 = Vm * Vm;
525: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
526: PD = QD = 0.0;
527: for (k = 0; k < ld_nsegsp[i]; k++) PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
528: for (k = 0; k < ld_nsegsq[i]; k++) QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
530: /* Load currents */
531: IDr = (PD * Vr + QD * Vi) / Vm2;
532: IDi = (-QD * Vr + PD * Vi) / Vm2;
534: fnet[2 * lbus[i]] += IDi;
535: fnet[2 * lbus[i] + 1] += IDr;
536: }
537: VecRestoreArray(user->V0, &v0);
539: VecRestoreArrayRead(Xgen, &xgen);
540: VecRestoreArrayRead(Xnet, &xnet);
541: VecRestoreArrayWrite(Fgen, &fgen);
542: VecRestoreArrayWrite(Fnet, &fnet);
544: DMCompositeGather(user->dmpgrid, INSERT_VALUES, F, Fgen, Fnet);
545: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
546: DMCompositeRestoreLocalVectors(user->dmpgrid, &Fgen, &Fnet);
547: return 0;
548: }
550: /* f(x,y)
551: g(x,y)
552: */
553: PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ctx)
554: {
555: Userctx *user = (Userctx *)ctx;
557: user->t = t;
558: ResidualFunction(X, F, user);
559: return 0;
560: }
562: /* f(x,y) - \dot{x}
563: g(x,y) = 0
564: */
565: PetscErrorCode IFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, void *ctx)
566: {
567: PetscScalar *f;
568: const PetscScalar *xdot;
569: PetscInt i;
571: RHSFunction(ts, t, X, F, ctx);
572: VecScale(F, -1.0);
573: VecGetArray(F, &f);
574: VecGetArrayRead(Xdot, &xdot);
575: for (i = 0; i < ngen; i++) {
576: f[9 * i] += xdot[9 * i];
577: f[9 * i + 1] += xdot[9 * i + 1];
578: f[9 * i + 2] += xdot[9 * i + 2];
579: f[9 * i + 3] += xdot[9 * i + 3];
580: f[9 * i + 6] += xdot[9 * i + 6];
581: f[9 * i + 7] += xdot[9 * i + 7];
582: f[9 * i + 8] += xdot[9 * i + 8];
583: }
584: VecRestoreArray(F, &f);
585: VecRestoreArrayRead(Xdot, &xdot);
586: return 0;
587: }
589: /* This function is used for solving the algebraic system only during fault on and
590: off times. It computes the entire F and then zeros out the part corresponding to
591: differential equations
592: F = [0;g(y)];
593: */
594: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
595: {
596: Userctx *user = (Userctx *)ctx;
597: PetscScalar *f;
598: PetscInt i;
600: ResidualFunction(X, F, user);
601: VecGetArray(F, &f);
602: for (i = 0; i < ngen; i++) {
603: f[9 * i] = 0;
604: f[9 * i + 1] = 0;
605: f[9 * i + 2] = 0;
606: f[9 * i + 3] = 0;
607: f[9 * i + 6] = 0;
608: f[9 * i + 7] = 0;
609: f[9 * i + 8] = 0;
610: }
611: VecRestoreArray(F, &f);
612: return 0;
613: }
615: PetscErrorCode PostStage(TS ts, PetscReal t, PetscInt i, Vec *X)
616: {
617: Userctx *user;
619: TSGetApplicationContext(ts, &user);
620: SNESSolve(user->snes_alg, NULL, X[i]);
621: return 0;
622: }
624: PetscErrorCode PostEvaluate(TS ts)
625: {
626: Userctx *user;
627: Vec X;
629: TSGetApplicationContext(ts, &user);
630: TSGetSolution(ts, &X);
631: SNESSolve(user->snes_alg, NULL, X);
632: return 0;
633: }
635: PetscErrorCode PreallocateJacobian(Mat J, Userctx *user)
636: {
637: PetscInt *d_nnz;
638: PetscInt i, idx = 0, start = 0;
639: PetscInt ncols;
641: PetscMalloc1(user->neqs_pgrid, &d_nnz);
642: for (i = 0; i < user->neqs_pgrid; i++) d_nnz[i] = 0;
643: /* Generator subsystem */
644: for (i = 0; i < ngen; i++) {
645: d_nnz[idx] += 3;
646: d_nnz[idx + 1] += 2;
647: d_nnz[idx + 2] += 2;
648: d_nnz[idx + 3] += 5;
649: d_nnz[idx + 4] += 6;
650: d_nnz[idx + 5] += 6;
652: d_nnz[user->neqs_gen + 2 * gbus[i]] += 3;
653: d_nnz[user->neqs_gen + 2 * gbus[i] + 1] += 3;
655: d_nnz[idx + 6] += 2;
656: d_nnz[idx + 7] += 2;
657: d_nnz[idx + 8] += 5;
659: idx = idx + 9;
660: }
661: start = user->neqs_gen;
663: for (i = 0; i < nbus; i++) {
664: MatGetRow(user->Ybus, 2 * i, &ncols, NULL, NULL);
665: d_nnz[start + 2 * i] += ncols;
666: d_nnz[start + 2 * i + 1] += ncols;
667: MatRestoreRow(user->Ybus, 2 * i, &ncols, NULL, NULL);
668: }
669: MatSeqAIJSetPreallocation(J, 0, d_nnz);
670: PetscFree(d_nnz);
671: return 0;
672: }
674: /*
675: J = [df_dx, df_dy
676: dg_dx, dg_dy]
677: */
678: PetscErrorCode ResidualJacobian(Vec X, Mat J, Mat B, void *ctx)
679: {
680: Userctx *user = (Userctx *)ctx;
681: Vec Xgen, Xnet;
682: const PetscScalar *xgen, *xnet;
683: PetscInt i, idx = 0;
684: PetscScalar Vr, Vi, Vm, Vm2;
685: PetscScalar Eqp, Edp, delta; /* Generator variables */
686: PetscScalar Efd;
687: PetscScalar Id, Iq; /* Generator dq axis currents */
688: PetscScalar Vd, Vq;
689: PetscScalar val[10];
690: PetscInt row[2], col[10];
691: PetscInt net_start = user->neqs_gen;
692: PetscScalar Zdq_inv[4], det;
693: PetscScalar dVd_dVr, dVd_dVi, dVq_dVr, dVq_dVi, dVd_ddelta, dVq_ddelta;
694: PetscScalar dIGr_ddelta, dIGi_ddelta, dIGr_dId, dIGr_dIq, dIGi_dId, dIGi_dIq;
695: PetscScalar dSE_dEfd;
696: PetscScalar dVm_dVd, dVm_dVq, dVm_dVr, dVm_dVi;
697: PetscInt ncols;
698: const PetscInt *cols;
699: const PetscScalar *yvals;
700: PetscInt k;
701: PetscScalar PD, QD, Vm0, Vm4;
702: const PetscScalar *v0;
703: PetscScalar dPD_dVr, dPD_dVi, dQD_dVr, dQD_dVi;
704: PetscScalar dIDr_dVr, dIDr_dVi, dIDi_dVr, dIDi_dVi;
706: MatZeroEntries(B);
707: DMCompositeGetLocalVectors(user->dmpgrid, &Xgen, &Xnet);
708: DMCompositeScatter(user->dmpgrid, X, Xgen, Xnet);
710: VecGetArrayRead(Xgen, &xgen);
711: VecGetArrayRead(Xnet, &xnet);
713: /* Generator subsystem */
714: for (i = 0; i < ngen; i++) {
715: Eqp = xgen[idx];
716: Edp = xgen[idx + 1];
717: delta = xgen[idx + 2];
718: Id = xgen[idx + 4];
719: Iq = xgen[idx + 5];
720: Efd = xgen[idx + 6];
722: /* fgen[idx] = (-Eqp - (Xd[i] - Xdp[i])*Id + Efd)/Td0p[i]; */
723: row[0] = idx;
724: col[0] = idx;
725: col[1] = idx + 4;
726: col[2] = idx + 6;
727: val[0] = -1 / Td0p[i];
728: val[1] = -(Xd[i] - Xdp[i]) / Td0p[i];
729: val[2] = 1 / Td0p[i];
731: MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);
733: /* fgen[idx+1] = (-Edp + (Xq[i] - Xqp[i])*Iq)/Tq0p[i]; */
734: row[0] = idx + 1;
735: col[0] = idx + 1;
736: col[1] = idx + 5;
737: val[0] = -1 / Tq0p[i];
738: val[1] = (Xq[i] - Xqp[i]) / Tq0p[i];
739: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
741: /* fgen[idx+2] = w - w_s; */
742: row[0] = idx + 2;
743: col[0] = idx + 2;
744: col[1] = idx + 3;
745: val[0] = 0;
746: val[1] = 1;
747: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
749: /* fgen[idx+3] = (TM[i] - Edp*Id - Eqp*Iq - (Xqp[i] - Xdp[i])*Id*Iq - D[i]*(w - w_s))/M[i]; */
750: row[0] = idx + 3;
751: col[0] = idx;
752: col[1] = idx + 1;
753: col[2] = idx + 3;
754: col[3] = idx + 4;
755: col[4] = idx + 5;
756: val[0] = -Iq / M[i];
757: val[1] = -Id / M[i];
758: val[2] = -D[i] / M[i];
759: val[3] = (-Edp - (Xqp[i] - Xdp[i]) * Iq) / M[i];
760: val[4] = (-Eqp - (Xqp[i] - Xdp[i]) * Id) / M[i];
761: MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES);
763: Vr = xnet[2 * gbus[i]]; /* Real part of generator terminal voltage */
764: Vi = xnet[2 * gbus[i] + 1]; /* Imaginary part of the generator terminal voltage */
765: ri2dq(Vr, Vi, delta, &Vd, &Vq);
767: det = Rs[i] * Rs[i] + Xdp[i] * Xqp[i];
769: Zdq_inv[0] = Rs[i] / det;
770: Zdq_inv[1] = Xqp[i] / det;
771: Zdq_inv[2] = -Xdp[i] / det;
772: Zdq_inv[3] = Rs[i] / det;
774: dVd_dVr = PetscSinScalar(delta);
775: dVd_dVi = -PetscCosScalar(delta);
776: dVq_dVr = PetscCosScalar(delta);
777: dVq_dVi = PetscSinScalar(delta);
778: dVd_ddelta = Vr * PetscCosScalar(delta) + Vi * PetscSinScalar(delta);
779: dVq_ddelta = -Vr * PetscSinScalar(delta) + Vi * PetscCosScalar(delta);
781: /* fgen[idx+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; */
782: row[0] = idx + 4;
783: col[0] = idx;
784: col[1] = idx + 1;
785: col[2] = idx + 2;
786: val[0] = -Zdq_inv[1];
787: val[1] = -Zdq_inv[0];
788: val[2] = Zdq_inv[0] * dVd_ddelta + Zdq_inv[1] * dVq_ddelta;
789: col[3] = idx + 4;
790: col[4] = net_start + 2 * gbus[i];
791: col[5] = net_start + 2 * gbus[i] + 1;
792: val[3] = 1;
793: val[4] = Zdq_inv[0] * dVd_dVr + Zdq_inv[1] * dVq_dVr;
794: val[5] = Zdq_inv[0] * dVd_dVi + Zdq_inv[1] * dVq_dVi;
795: MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES);
797: /* fgen[idx+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; */
798: row[0] = idx + 5;
799: col[0] = idx;
800: col[1] = idx + 1;
801: col[2] = idx + 2;
802: val[0] = -Zdq_inv[3];
803: val[1] = -Zdq_inv[2];
804: val[2] = Zdq_inv[2] * dVd_ddelta + Zdq_inv[3] * dVq_ddelta;
805: col[3] = idx + 5;
806: col[4] = net_start + 2 * gbus[i];
807: col[5] = net_start + 2 * gbus[i] + 1;
808: val[3] = 1;
809: val[4] = Zdq_inv[2] * dVd_dVr + Zdq_inv[3] * dVq_dVr;
810: val[5] = Zdq_inv[2] * dVd_dVi + Zdq_inv[3] * dVq_dVi;
811: MatSetValues(J, 1, row, 6, col, val, INSERT_VALUES);
813: dIGr_ddelta = Id * PetscCosScalar(delta) - Iq * PetscSinScalar(delta);
814: dIGi_ddelta = Id * PetscSinScalar(delta) + Iq * PetscCosScalar(delta);
815: dIGr_dId = PetscSinScalar(delta);
816: dIGr_dIq = PetscCosScalar(delta);
817: dIGi_dId = -PetscCosScalar(delta);
818: dIGi_dIq = PetscSinScalar(delta);
820: /* fnet[2*gbus[i]] -= IGi; */
821: row[0] = net_start + 2 * gbus[i];
822: col[0] = idx + 2;
823: col[1] = idx + 4;
824: col[2] = idx + 5;
825: val[0] = -dIGi_ddelta;
826: val[1] = -dIGi_dId;
827: val[2] = -dIGi_dIq;
828: MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);
830: /* fnet[2*gbus[i]+1] -= IGr; */
831: row[0] = net_start + 2 * gbus[i] + 1;
832: col[0] = idx + 2;
833: col[1] = idx + 4;
834: col[2] = idx + 5;
835: val[0] = -dIGr_ddelta;
836: val[1] = -dIGr_dId;
837: val[2] = -dIGr_dIq;
838: MatSetValues(J, 1, row, 3, col, val, INSERT_VALUES);
840: Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
842: /* fgen[idx+6] = (-KE[i]*Efd - SE + VR)/TE[i]; */
843: /* SE = k1[i]*PetscExpScalar(k2[i]*Efd); */
845: dSE_dEfd = k1[i] * k2[i] * PetscExpScalar(k2[i] * Efd);
847: row[0] = idx + 6;
848: col[0] = idx + 6;
849: col[1] = idx + 8;
850: val[0] = (-KE[i] - dSE_dEfd) / TE[i];
851: val[1] = 1 / TE[i];
852: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
854: /* Exciter differential equations */
856: /* fgen[idx+7] = (-RF + KF[i]*Efd/TF[i])/TF[i]; */
857: row[0] = idx + 7;
858: col[0] = idx + 6;
859: col[1] = idx + 7;
860: val[0] = (KF[i] / TF[i]) / TF[i];
861: val[1] = -1 / TF[i];
862: MatSetValues(J, 1, row, 2, col, val, INSERT_VALUES);
864: /* fgen[idx+8] = (-VR + KA[i]*RF - KA[i]*KF[i]*Efd/TF[i] + KA[i]*(Vref[i] - Vm))/TA[i]; */
865: /* Vm = (Vd^2 + Vq^2)^0.5; */
867: row[0] = idx + 8;
868: if (VRatmax[i]) {
869: col[0] = idx + 8;
870: val[0] = 1.0;
871: MatSetValues(J, 1, row, 1, col, val, INSERT_VALUES);
872: } else if (VRatmin[i]) {
873: col[0] = idx + 8;
874: val[0] = -1.0;
875: MatSetValues(J, 1, row, 1, col, val, INSERT_VALUES);
876: } else {
877: dVm_dVd = Vd / Vm;
878: dVm_dVq = Vq / Vm;
879: dVm_dVr = dVm_dVd * dVd_dVr + dVm_dVq * dVq_dVr;
880: dVm_dVi = dVm_dVd * dVd_dVi + dVm_dVq * dVq_dVi;
881: row[0] = idx + 8;
882: col[0] = idx + 6;
883: col[1] = idx + 7;
884: col[2] = idx + 8;
885: val[0] = -(KA[i] * KF[i] / TF[i]) / TA[i];
886: val[1] = KA[i] / TA[i];
887: val[2] = -1 / TA[i];
888: col[3] = net_start + 2 * gbus[i];
889: col[4] = net_start + 2 * gbus[i] + 1;
890: val[3] = -KA[i] * dVm_dVr / TA[i];
891: val[4] = -KA[i] * dVm_dVi / TA[i];
892: MatSetValues(J, 1, row, 5, col, val, INSERT_VALUES);
893: }
894: idx = idx + 9;
895: }
897: for (i = 0; i < nbus; i++) {
898: MatGetRow(user->Ybus, 2 * i, &ncols, &cols, &yvals);
899: row[0] = net_start + 2 * i;
900: for (k = 0; k < ncols; k++) {
901: col[k] = net_start + cols[k];
902: val[k] = yvals[k];
903: }
904: MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES);
905: MatRestoreRow(user->Ybus, 2 * i, &ncols, &cols, &yvals);
907: MatGetRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals);
908: row[0] = net_start + 2 * i + 1;
909: for (k = 0; k < ncols; k++) {
910: col[k] = net_start + cols[k];
911: val[k] = yvals[k];
912: }
913: MatSetValues(J, 1, row, ncols, col, val, INSERT_VALUES);
914: MatRestoreRow(user->Ybus, 2 * i + 1, &ncols, &cols, &yvals);
915: }
917: MatAssemblyBegin(J, MAT_FLUSH_ASSEMBLY);
918: MatAssemblyEnd(J, MAT_FLUSH_ASSEMBLY);
920: VecGetArrayRead(user->V0, &v0);
921: for (i = 0; i < nload; i++) {
922: Vr = xnet[2 * lbus[i]]; /* Real part of load bus voltage */
923: Vi = xnet[2 * lbus[i] + 1]; /* Imaginary part of the load bus voltage */
924: Vm = PetscSqrtScalar(Vr * Vr + Vi * Vi);
925: Vm2 = Vm * Vm;
926: Vm4 = Vm2 * Vm2;
927: Vm0 = PetscSqrtScalar(v0[2 * lbus[i]] * v0[2 * lbus[i]] + v0[2 * lbus[i] + 1] * v0[2 * lbus[i] + 1]);
928: PD = QD = 0.0;
929: dPD_dVr = dPD_dVi = dQD_dVr = dQD_dVi = 0.0;
930: for (k = 0; k < ld_nsegsp[i]; k++) {
931: PD += ld_alphap[k] * PD0[i] * PetscPowScalar((Vm / Vm0), ld_betap[k]);
932: dPD_dVr += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vr * PetscPowScalar(Vm, (ld_betap[k] - 2));
933: dPD_dVi += ld_alphap[k] * ld_betap[k] * PD0[i] * PetscPowScalar((1 / Vm0), ld_betap[k]) * Vi * PetscPowScalar(Vm, (ld_betap[k] - 2));
934: }
935: for (k = 0; k < ld_nsegsq[i]; k++) {
936: QD += ld_alphaq[k] * QD0[i] * PetscPowScalar((Vm / Vm0), ld_betaq[k]);
937: dQD_dVr += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vr * PetscPowScalar(Vm, (ld_betaq[k] - 2));
938: dQD_dVi += ld_alphaq[k] * ld_betaq[k] * QD0[i] * PetscPowScalar((1 / Vm0), ld_betaq[k]) * Vi * PetscPowScalar(Vm, (ld_betaq[k] - 2));
939: }
941: /* IDr = (PD*Vr + QD*Vi)/Vm2; */
942: /* IDi = (-QD*Vr + PD*Vi)/Vm2; */
944: dIDr_dVr = (dPD_dVr * Vr + dQD_dVr * Vi + PD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vr) / Vm4;
945: dIDr_dVi = (dPD_dVi * Vr + dQD_dVi * Vi + QD) / Vm2 - ((PD * Vr + QD * Vi) * 2 * Vi) / Vm4;
947: dIDi_dVr = (-dQD_dVr * Vr + dPD_dVr * Vi - QD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vr) / Vm4;
948: dIDi_dVi = (-dQD_dVi * Vr + dPD_dVi * Vi + PD) / Vm2 - ((-QD * Vr + PD * Vi) * 2 * Vi) / Vm4;
950: /* fnet[2*lbus[i]] += IDi; */
951: row[0] = net_start + 2 * lbus[i];
952: col[0] = net_start + 2 * lbus[i];
953: col[1] = net_start + 2 * lbus[i] + 1;
954: val[0] = dIDi_dVr;
955: val[1] = dIDi_dVi;
956: MatSetValues(J, 1, row, 2, col, val, ADD_VALUES);
957: /* fnet[2*lbus[i]+1] += IDr; */
958: row[0] = net_start + 2 * lbus[i] + 1;
959: col[0] = net_start + 2 * lbus[i];
960: col[1] = net_start + 2 * lbus[i] + 1;
961: val[0] = dIDr_dVr;
962: val[1] = dIDr_dVi;
963: MatSetValues(J, 1, row, 2, col, val, ADD_VALUES);
964: }
965: VecRestoreArrayRead(user->V0, &v0);
967: VecRestoreArrayRead(Xgen, &xgen);
968: VecRestoreArrayRead(Xnet, &xnet);
970: DMCompositeRestoreLocalVectors(user->dmpgrid, &Xgen, &Xnet);
972: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
973: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
974: return 0;
975: }
977: /*
978: J = [I, 0
979: dg_dx, dg_dy]
980: */
981: PetscErrorCode AlgJacobian(SNES snes, Vec X, Mat A, Mat B, void *ctx)
982: {
983: Userctx *user = (Userctx *)ctx;
985: ResidualJacobian(X, A, B, ctx);
986: MatSetOption(A, MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE);
987: MatZeroRowsIS(A, user->is_diff, 1.0, NULL, NULL);
988: return 0;
989: }
991: /*
992: J = [-df_dx, -df_dy
993: dg_dx, dg_dy]
994: */
996: PetscErrorCode RHSJacobian(TS ts, PetscReal t, Vec X, Mat A, Mat B, void *ctx)
997: {
998: Userctx *user = (Userctx *)ctx;
1000: user->t = t;
1002: ResidualJacobian(X, A, B, user);
1004: return 0;
1005: }
1007: /*
1008: J = [df_dx-aI, df_dy
1009: dg_dx, dg_dy]
1010: */
1012: PetscErrorCode IJacobian(TS ts, PetscReal t, Vec X, Vec Xdot, PetscReal a, Mat A, Mat B, Userctx *user)
1013: {
1014: PetscScalar atmp = (PetscScalar)a;
1015: PetscInt i, row;
1017: user->t = t;
1019: RHSJacobian(ts, t, X, A, B, user);
1020: MatScale(B, -1.0);
1021: for (i = 0; i < ngen; i++) {
1022: row = 9 * i;
1023: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
1024: row = 9 * i + 1;
1025: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
1026: row = 9 * i + 2;
1027: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
1028: row = 9 * i + 3;
1029: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
1030: row = 9 * i + 6;
1031: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
1032: row = 9 * i + 7;
1033: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
1034: row = 9 * i + 8;
1035: MatSetValues(A, 1, &row, 1, &row, &atmp, ADD_VALUES);
1036: }
1037: MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
1038: MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);
1039: return 0;
1040: }
1042: int main(int argc, char **argv)
1043: {
1044: TS ts;
1045: SNES snes_alg;
1046: PetscMPIInt size;
1047: Userctx user;
1048: PetscViewer Xview, Ybusview, viewer;
1049: Vec X, F_alg;
1050: Mat J, A;
1051: PetscInt i, idx, *idx2;
1052: Vec Xdot;
1053: PetscScalar *x, *mat, *amat;
1054: const PetscScalar *rmat;
1055: Vec vatol;
1056: PetscInt *direction;
1057: PetscBool *terminate;
1058: const PetscInt *idx3;
1059: PetscScalar *vatoli;
1060: PetscInt k;
1063: PetscInitialize(&argc, &argv, "petscoptions", help);
1064: MPI_Comm_size(PETSC_COMM_WORLD, &size);
1067: user.neqs_gen = 9 * ngen; /* # eqs. for generator subsystem */
1068: user.neqs_net = 2 * nbus; /* # eqs. for network subsystem */
1069: user.neqs_pgrid = user.neqs_gen + user.neqs_net;
1071: /* Create indices for differential and algebraic equations */
1073: PetscMalloc1(7 * ngen, &idx2);
1074: for (i = 0; i < ngen; i++) {
1075: idx2[7 * i] = 9 * i;
1076: idx2[7 * i + 1] = 9 * i + 1;
1077: idx2[7 * i + 2] = 9 * i + 2;
1078: idx2[7 * i + 3] = 9 * i + 3;
1079: idx2[7 * i + 4] = 9 * i + 6;
1080: idx2[7 * i + 5] = 9 * i + 7;
1081: idx2[7 * i + 6] = 9 * i + 8;
1082: }
1083: ISCreateGeneral(PETSC_COMM_WORLD, 7 * ngen, idx2, PETSC_COPY_VALUES, &user.is_diff);
1084: ISComplement(user.is_diff, 0, user.neqs_pgrid, &user.is_alg);
1085: PetscFree(idx2);
1087: /* Read initial voltage vector and Ybus */
1088: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "X.bin", FILE_MODE_READ, &Xview);
1089: PetscViewerBinaryOpen(PETSC_COMM_WORLD, "Ybus.bin", FILE_MODE_READ, &Ybusview);
1091: VecCreate(PETSC_COMM_WORLD, &user.V0);
1092: VecSetSizes(user.V0, PETSC_DECIDE, user.neqs_net);
1093: VecLoad(user.V0, Xview);
1095: MatCreate(PETSC_COMM_WORLD, &user.Ybus);
1096: MatSetSizes(user.Ybus, PETSC_DECIDE, PETSC_DECIDE, user.neqs_net, user.neqs_net);
1097: MatSetType(user.Ybus, MATBAIJ);
1098: /* MatSetBlockSize(user.Ybus,2); */
1099: MatLoad(user.Ybus, Ybusview);
1101: /* Set run time options */
1102: PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1103: {
1104: user.tfaulton = 1.0;
1105: user.tfaultoff = 1.2;
1106: user.Rfault = 0.0001;
1107: user.setisdiff = PETSC_FALSE;
1108: user.semiexplicit = PETSC_FALSE;
1109: user.faultbus = 8;
1110: PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL);
1111: PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL);
1112: PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL);
1113: user.t0 = 0.0;
1114: user.tmax = 5.0;
1115: PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL);
1116: PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL);
1117: PetscOptionsBool("-setisdiff", "", "", user.setisdiff, &user.setisdiff, NULL);
1118: PetscOptionsBool("-dae_semiexplicit", "", "", user.semiexplicit, &user.semiexplicit, NULL);
1119: }
1120: PetscOptionsEnd();
1122: PetscViewerDestroy(&Xview);
1123: PetscViewerDestroy(&Ybusview);
1125: /* Create DMs for generator and network subsystems */
1126: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_gen, 1, 1, NULL, &user.dmgen);
1127: DMSetOptionsPrefix(user.dmgen, "dmgen_");
1128: DMSetFromOptions(user.dmgen);
1129: DMSetUp(user.dmgen);
1130: DMDACreate1d(PETSC_COMM_WORLD, DM_BOUNDARY_NONE, user.neqs_net, 1, 1, NULL, &user.dmnet);
1131: DMSetOptionsPrefix(user.dmnet, "dmnet_");
1132: DMSetFromOptions(user.dmnet);
1133: DMSetUp(user.dmnet);
1134: /* Create a composite DM packer and add the two DMs */
1135: DMCompositeCreate(PETSC_COMM_WORLD, &user.dmpgrid);
1136: DMSetOptionsPrefix(user.dmpgrid, "pgrid_");
1137: DMCompositeAddDM(user.dmpgrid, user.dmgen);
1138: DMCompositeAddDM(user.dmpgrid, user.dmnet);
1140: DMCreateGlobalVector(user.dmpgrid, &X);
1142: MatCreate(PETSC_COMM_WORLD, &J);
1143: MatSetSizes(J, PETSC_DECIDE, PETSC_DECIDE, user.neqs_pgrid, user.neqs_pgrid);
1144: MatSetFromOptions(J);
1145: PreallocateJacobian(J, &user);
1147: /* Create matrix to save solutions at each time step */
1148: user.stepnum = 0;
1150: MatCreateSeqDense(PETSC_COMM_SELF, user.neqs_pgrid + 1, 1002, NULL, &user.Sol);
1151: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1152: Create timestepping solver context
1153: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1154: TSCreate(PETSC_COMM_WORLD, &ts);
1155: TSSetProblemType(ts, TS_NONLINEAR);
1156: if (user.semiexplicit) {
1157: TSSetType(ts, TSRK);
1158: TSSetRHSFunction(ts, NULL, RHSFunction, &user);
1159: TSSetRHSJacobian(ts, J, J, RHSJacobian, &user);
1160: } else {
1161: TSSetType(ts, TSCN);
1162: TSSetEquationType(ts, TS_EQ_DAE_IMPLICIT_INDEX1);
1163: TSSetIFunction(ts, NULL, (TSIFunction)IFunction, &user);
1164: TSSetIJacobian(ts, J, J, (TSIJacobian)IJacobian, &user);
1165: }
1166: TSSetApplicationContext(ts, &user);
1168: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
1169: Set initial conditions
1170: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
1171: SetInitialGuess(X, &user);
1172: /* Just to set up the Jacobian structure */
1174: VecDuplicate(X, &Xdot);
1175: IJacobian(ts, 0.0, X, Xdot, 0.0, J, J, &user);
1176: VecDestroy(&Xdot);
1178: /* Save initial solution */
1180: idx = user.stepnum * (user.neqs_pgrid + 1);
1181: MatDenseGetArray(user.Sol, &mat);
1182: VecGetArray(X, &x);
1184: mat[idx] = 0.0;
1186: PetscArraycpy(mat + idx + 1, x, user.neqs_pgrid);
1187: MatDenseRestoreArray(user.Sol, &mat);
1188: VecRestoreArray(X, &x);
1189: user.stepnum++;
1191: TSSetMaxTime(ts, user.tmax);
1192: TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP);
1193: TSSetTimeStep(ts, 0.01);
1194: TSSetFromOptions(ts);
1195: TSSetPostStep(ts, SaveSolution);
1196: TSSetSolution(ts, X);
1198: PetscMalloc1((2 * ngen + 2), &direction);
1199: PetscMalloc1((2 * ngen + 2), &terminate);
1200: direction[0] = direction[1] = 1;
1201: terminate[0] = terminate[1] = PETSC_FALSE;
1202: for (i = 0; i < ngen; i++) {
1203: direction[2 + 2 * i] = -1;
1204: direction[2 + 2 * i + 1] = 1;
1205: terminate[2 + 2 * i] = terminate[2 + 2 * i + 1] = PETSC_FALSE;
1206: }
1208: TSSetEventHandler(ts, 2 * ngen + 2, direction, terminate, EventFunction, PostEventFunction, (void *)&user);
1210: if (user.semiexplicit) {
1211: /* Use a semi-explicit approach with the time-stepping done by an explicit method and the
1212: algrebraic part solved via PostStage and PostEvaluate callbacks
1213: */
1214: TSSetType(ts, TSRK);
1215: TSSetPostStage(ts, PostStage);
1216: TSSetPostEvaluate(ts, PostEvaluate);
1217: }
1219: if (user.setisdiff) {
1220: /* Create vector of absolute tolerances and set the algebraic part to infinity */
1221: VecDuplicate(X, &vatol);
1222: VecSet(vatol, 100000.0);
1223: VecGetArray(vatol, &vatoli);
1224: ISGetIndices(user.is_diff, &idx3);
1225: for (k = 0; k < 7 * ngen; k++) vatoli[idx3[k]] = 1e-2;
1226: VecRestoreArray(vatol, &vatoli);
1227: }
1229: /* Create the nonlinear solver for solving the algebraic system */
1230: /* Note that although the algebraic system needs to be solved only for
1231: Idq and V, we reuse the entire system including xgen. The xgen
1232: variables are held constant by setting their residuals to 0 and
1233: putting a 1 on the Jacobian diagonal for xgen rows
1234: */
1236: VecDuplicate(X, &F_alg);
1237: SNESCreate(PETSC_COMM_WORLD, &snes_alg);
1238: SNESSetFunction(snes_alg, F_alg, AlgFunction, &user);
1239: SNESSetJacobian(snes_alg, J, J, AlgJacobian, &user);
1240: SNESSetFromOptions(snes_alg);
1242: user.snes_alg = snes_alg;
1244: /* Solve */
1245: TSSolve(ts, X);
1247: MatAssemblyBegin(user.Sol, MAT_FINAL_ASSEMBLY);
1248: MatAssemblyEnd(user.Sol, MAT_FINAL_ASSEMBLY);
1250: MatCreateSeqDense(PETSC_COMM_SELF, user.neqs_pgrid + 1, user.stepnum, NULL, &A);
1251: MatDenseGetArrayRead(user.Sol, &rmat);
1252: MatDenseGetArray(A, &amat);
1253: PetscArraycpy(amat, rmat, user.stepnum * (user.neqs_pgrid + 1));
1254: MatDenseRestoreArray(A, &amat);
1255: MatDenseRestoreArrayRead(user.Sol, &rmat);
1256: PetscViewerBinaryOpen(PETSC_COMM_SELF, "out.bin", FILE_MODE_WRITE, &viewer);
1257: MatView(A, viewer);
1258: PetscViewerDestroy(&viewer);
1259: MatDestroy(&A);
1261: PetscFree(direction);
1262: PetscFree(terminate);
1263: SNESDestroy(&snes_alg);
1264: VecDestroy(&F_alg);
1265: MatDestroy(&J);
1266: MatDestroy(&user.Ybus);
1267: MatDestroy(&user.Sol);
1268: VecDestroy(&X);
1269: VecDestroy(&user.V0);
1270: DMDestroy(&user.dmgen);
1271: DMDestroy(&user.dmnet);
1272: DMDestroy(&user.dmpgrid);
1273: ISDestroy(&user.is_diff);
1274: ISDestroy(&user.is_alg);
1275: TSDestroy(&ts);
1276: if (user.setisdiff) VecDestroy(&vatol);
1277: PetscFinalize();
1278: return 0;
1279: }
1281: /*TEST
1283: build:
1284: requires: double !complex !defined(PETSC_USE_64BIT_INDICES)
1286: test:
1287: suffix: implicit
1288: args: -ts_monitor -snes_monitor_short
1289: localrunfiles: petscoptions X.bin Ybus.bin
1291: test:
1292: suffix: semiexplicit
1293: args: -ts_monitor -snes_monitor_short -dae_semiexplicit -ts_rk_type 2a
1294: localrunfiles: petscoptions X.bin Ybus.bin
1296: test:
1297: suffix: steprestart
1298: # needs ARKIMEX methods with all implicit stages since the mass matrix is not the identity
1299: args: -ts_monitor -snes_monitor_short -ts_type arkimex -ts_arkimex_type prssp2
1300: localrunfiles: petscoptions X.bin Ybus.bin
1302: TEST*/