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;

158:   PetscFunctionBeginUser;
159:   /* Read V0 and Ybus from files */
160:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "X.bin", FILE_MODE_READ, &Xview));
161:   PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF, "Ybus.bin", FILE_MODE_READ, &Ybusview));
162:   PetscCall(VecCreate(PETSC_COMM_SELF, &V0));
163:   PetscCall(VecLoad(V0, Xview));

165:   PetscCall(MatCreate(PETSC_COMM_SELF, &Ybus));
166:   PetscCall(MatSetType(Ybus, MATBAIJ));
167:   PetscCall(MatLoad(Ybus, Ybusview));

169:   /* Destroy unnecessary stuff */
170:   PetscCall(PetscViewerDestroy(&Xview));
171:   PetscCall(PetscViewerDestroy(&Ybusview));

173:   PetscCall(MatGetLocalSize(Ybus, &m, &n));
174:   neqs_net = 2 * NBUS; /* # eqs. for network subsystem   */
175:   PetscCheck(m == neqs_net && n == neqs_net, PETSC_COMM_SELF, PETSC_ERR_ARG_SIZ, "matrix Ybus is in wrong sizes");

177:   M[0] = 2 * H[0] / W_S;
178:   M[1] = 2 * H[1] / W_S;
179:   M[2] = 2 * H[2] / W_S;
180:   D[0] = 0.1 * M[0];
181:   D[1] = 0.1 * M[1];
182:   D[2] = 0.1 * M[2];

184:   /* Allocate memory for bus, generators, exciter, loads and branches */
185:   PetscCall(PetscCalloc5(NBUS * nc, &bus, NGEN * nc, &gen, NLOAD * nc, &load, NBRANCH * nc + (nc - 1), &branch, NGEN * nc, &exc));

187:   PetscCall(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:       PetscCall(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:   PetscCall(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:       PetscCall(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:   PetscCall(VecRestoreArrayRead(V0, &varr));

351:   /* Destroy unnecessary stuff */
352:   PetscCall(MatDestroy(&Ybus));
353:   PetscCall(VecDestroy(&V0));
354:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
378:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));
379:   PetscCall(DMGetLocalVector(networkdm, &localX));

381:   PetscCall(VecSet(X, 0.0));
382:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
383:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

385:   PetscCall(VecGetArray(localX, &xarr));

387:   for (v = vStart; v < vEnd; v++) {
388:     PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
389:     if (ghostvtex) continue;

391:     PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));
392:     for (j = 0; j < numComps; j++) {
393:       PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
394:       if (key == 1) {
395:         bus = (Bus *)(component);

397:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
398:         xarr[offset]     = bus->vr;
399:         xarr[offset + 1] = bus->vi;

401:         Vr = bus->vr;
402:         Vi = bus->vi;
403:       } else if (key == 2) {
404:         gen = (Gen *)(component);
405:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));
406:         Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
407:         Vm2 = Vm * Vm;
408:         /* Real part of gen current */
409:         IGr = (Vr * gen->PG + Vi * gen->QG) / Vm2;
410:         /* Imaginary part of gen current */
411:         IGi = (Vi * gen->PG - Vr * gen->QG) / Vm2;

413:         /* Machine angle */
414:         delta = atan2(Vi + gen->Xq * IGr, Vr - gen->Xq * IGi);
415:         theta = PETSC_PI / 2.0 - delta;

417:         /* d-axis stator current */
418:         Id = IGr * PetscCosScalar(theta) - IGi * PetscSinScalar(theta);

420:         /* q-axis stator current */
421:         Iq = IGr * PetscSinScalar(theta) + IGi * PetscCosScalar(theta);

423:         Vd = Vr * PetscCosScalar(theta) - Vi * PetscSinScalar(theta);
424:         Vq = Vr * PetscSinScalar(theta) + Vi * PetscCosScalar(theta);

426:         /* d-axis transient EMF */
427:         Edp = Vd + gen->Rs * Id - gen->Xqp * Iq;

429:         /* q-axis transient EMF */
430:         Eqp = Vq + gen->Rs * Iq + gen->Xdp * Id;

432:         gen->TM = gen->PG;

434:         xarr[offset]     = Eqp;
435:         xarr[offset + 1] = Edp;
436:         xarr[offset + 2] = delta;
437:         xarr[offset + 3] = W_S;
438:         xarr[offset + 4] = Id;
439:         xarr[offset + 5] = Iq;

441:         Efd = Eqp + (gen->Xd - gen->Xdp) * Id;

443:       } else if (key == 3) {
444:         exc = (Exc *)(component);
445:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offset));

447:         SE = exc->k1 * PetscExpScalar(exc->k2 * Efd);
448:         VR = exc->KE * Efd + SE;
449:         RF = exc->KF * Efd / exc->TF;

451:         xarr[offset]     = Efd;
452:         xarr[offset + 1] = RF;
453:         xarr[offset + 2] = VR;

455:         exc->Vref = Vm + (VR / exc->KA);
456:       }
457:     }
458:   }
459:   PetscCall(VecRestoreArray(localX, &xarr));
460:   PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
461:   PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
462:   PetscCall(DMRestoreLocalVector(networkdm, &localX));
463:   PetscFunctionReturn(PETSC_SUCCESS);
464: }

466: /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */
467: PetscErrorCode dq2ri(PetscScalar Fd, PetscScalar Fq, PetscScalar delta, PetscScalar *Fr, PetscScalar *Fi)
468: {
469:   PetscFunctionBegin;
470:   *Fr = Fd * PetscSinScalar(delta) + Fq * PetscCosScalar(delta);
471:   *Fi = -Fd * PetscCosScalar(delta) + Fq * PetscSinScalar(delta);
472:   PetscFunctionReturn(PETSC_SUCCESS);
473: }

475: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
476: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
477: {
478:   PetscFunctionBegin;
479:   *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
480:   *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
481:   PetscFunctionReturn(PETSC_SUCCESS);
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:   PetscFunctionBegin;
498:   PetscCall(VecSet(F, 0.0));

500:   PetscCall(TSGetDM(ts, &networkdm));
501:   PetscCall(DMGetLocalVector(networkdm, &localF));
502:   PetscCall(DMGetLocalVector(networkdm, &localX));
503:   PetscCall(DMGetLocalVector(networkdm, &localXdot));
504:   PetscCall(VecSet(localF, 0.0));

506:   /* update ghost values of localX and localXdot */
507:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
508:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

510:   PetscCall(DMGlobalToLocalBegin(networkdm, Xdot, INSERT_VALUES, localXdot));
511:   PetscCall(DMGlobalToLocalEnd(networkdm, Xdot, INSERT_VALUES, localXdot));

513:   PetscCall(VecGetArrayRead(localX, &xarr));
514:   PetscCall(VecGetArrayRead(localXdot, &xdotarr));
515:   PetscCall(VecGetArray(localF, &farr));

517:   PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

519:   for (v = vStart; v < vEnd; v++) {
520:     PetscInt    i, j, offsetbus, offsetgen, offsetexc, key;
521:     Bus        *bus;
522:     Gen        *gen;
523:     Exc        *exc;
524:     Load       *load;
525:     PetscBool   ghostvtex;
526:     PetscInt    numComps;
527:     PetscScalar Yffr, Yffi; /* Real and imaginary fault admittances */
528:     PetscScalar Vm, Vm2, Vm0;
529:     PetscScalar Vr0 = 0, Vi0 = 0;
530:     PetscScalar PD, QD;

532:     PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
533:     PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));

535:     for (j = 0; j < numComps; j++) {
536:       PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
537:       if (key == 1) {
538:         PetscInt        nconnedges;
539:         const PetscInt *connedges;

541:         bus = (Bus *)(component);
542:         PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus));
543:         if (!ghostvtex) {
544:           Vr = xarr[offsetbus];
545:           Vi = xarr[offsetbus + 1];

547:           Yffr = bus->yff[1];
548:           Yffi = bus->yff[0];

550:           if (user->alg_flg) {
551:             Yffr += user->ybusfault[bus->id * 2 + 1];
552:             Yffi += user->ybusfault[bus->id * 2];
553:           }
554:           Vr0 = bus->vr;
555:           Vi0 = bus->vi;

557:           /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
558:            The generator current injection, IG, and load current injection, ID are added later
559:            */
560:           farr[offsetbus] += Yffi * Vr + Yffr * Vi;     /* imaginary current due to diagonal elements */
561:           farr[offsetbus + 1] += Yffr * Vr - Yffi * Vi; /* real current due to diagonal elements */
562:         }

564:         PetscCall(DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges));

566:         for (i = 0; i < nconnedges; i++) {
567:           Branch         *branch;
568:           PetscInt        keye;
569:           PetscScalar     Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
570:           const PetscInt *cone;

572:           e = connedges[i];
573:           PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));

575:           Yfti = branch->yft[0];
576:           Yftr = branch->yft[1];

578:           PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));

580:           vfrom = cone[0];
581:           vto   = cone[1];

583:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
584:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));

586:           /* From bus and to bus real and imaginary voltages */
587:           Vfr = xarr[offsetfrom];
588:           Vfi = xarr[offsetfrom + 1];
589:           Vtr = xarr[offsetto];
590:           Vti = xarr[offsetto + 1];

592:           if (vfrom == v) {
593:             farr[offsetfrom] += Yftr * Vti + Yfti * Vtr;
594:             farr[offsetfrom + 1] += Yftr * Vtr - Yfti * Vti;
595:           } else {
596:             farr[offsetto] += Yftr * Vfi + Yfti * Vfr;
597:             farr[offsetto + 1] += Yftr * Vfr - Yfti * Vfi;
598:           }
599:         }
600:       } else if (key == 2) {
601:         if (!ghostvtex) {
602:           PetscScalar Eqp, Edp, delta, w; /* Generator variables */
603:           PetscScalar Efd;                /* Exciter field voltage */
604:           PetscScalar Id, Iq;             /* Generator dq axis currents */
605:           PetscScalar IGr, IGi, Zdq_inv[4], det;
606:           PetscScalar Xd, Xdp, Td0p, Xq, Xqp, Tq0p, TM, D, M, Rs; /* Generator parameters */

608:           gen = (Gen *)(component);
609:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));

611:           /* Generator state variables */
612:           Eqp   = xarr[offsetgen];
613:           Edp   = xarr[offsetgen + 1];
614:           delta = xarr[offsetgen + 2];
615:           w     = xarr[offsetgen + 3];
616:           Id    = xarr[offsetgen + 4];
617:           Iq    = xarr[offsetgen + 5];

619:           /* Generator parameters */
620:           Xd   = gen->Xd;
621:           Xdp  = gen->Xdp;
622:           Td0p = gen->Td0p;
623:           Xq   = gen->Xq;
624:           Xqp  = gen->Xqp;
625:           Tq0p = gen->Tq0p;
626:           TM   = gen->TM;
627:           D    = gen->D;
628:           M    = gen->M;
629:           Rs   = gen->Rs;

631:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, 2, &offsetexc));
632:           Efd = xarr[offsetexc];

634:           /* Generator differential equations */
635:           farr[offsetgen]     = (Eqp + (Xd - Xdp) * Id - Efd) / Td0p + xdotarr[offsetgen];
636:           farr[offsetgen + 1] = (Edp - (Xq - Xqp) * Iq) / Tq0p + xdotarr[offsetgen + 1];
637:           farr[offsetgen + 2] = -w + W_S + xdotarr[offsetgen + 2];
638:           farr[offsetgen + 3] = (-TM + Edp * Id + Eqp * Iq + (Xqp - Xdp) * Id * Iq + D * (w - W_S)) / M + xdotarr[offsetgen + 3];

640:           PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));

642:           /* Algebraic equations for stator currents */
643:           det = Rs * Rs + Xdp * Xqp;

645:           Zdq_inv[0] = Rs / det;
646:           Zdq_inv[1] = Xqp / det;
647:           Zdq_inv[2] = -Xdp / det;
648:           Zdq_inv[3] = Rs / det;

650:           farr[offsetgen + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
651:           farr[offsetgen + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;

653:           PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));

655:           /* Add generator current injection to network */
656:           farr[offsetbus] -= IGi;
657:           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:           PetscCall(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];
688:         }
689:       } else if (key == 4) {
690:         if (!ghostvtex) {
691:           PetscInt     k;
692:           PetscInt     ld_nsegsp;
693:           PetscInt     ld_nsegsq;
694:           PetscScalar *ld_alphap;
695:           PetscScalar *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;

697:           load = (Load *)(component);

699:           /* Load Parameters */
700:           ld_nsegsp = load->ld_nsegsp;
701:           ld_alphap = load->ld_alphap;
702:           ld_betap  = load->ld_betap;
703:           ld_nsegsq = load->ld_nsegsq;
704:           ld_alphaq = load->ld_alphaq;
705:           ld_betaq  = load->ld_betaq;
706:           PD0       = load->PD0;
707:           QD0       = load->QD0;

709:           Vr  = xarr[offsetbus];     /* Real part of generator terminal voltage */
710:           Vi  = xarr[offsetbus + 1]; /* Imaginary part of the generator terminal voltage */
711:           Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
712:           Vm2 = Vm * Vm;
713:           Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
714:           PD = QD = 0.0;
715:           for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar((Vm / Vm0), ld_betap[k]);
716:           for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar((Vm / Vm0), ld_betaq[k]);

718:           /* Load currents */
719:           IDr = (PD * Vr + QD * Vi) / Vm2;
720:           IDi = (-QD * Vr + PD * Vi) / Vm2;

722:           /* Load current contribution to the network */
723:           farr[offsetbus] += IDi;
724:           farr[offsetbus + 1] += IDr;
725:         }
726:       }
727:     }
728:   }

730:   PetscCall(VecRestoreArrayRead(localX, &xarr));
731:   PetscCall(VecRestoreArrayRead(localXdot, &xdotarr));
732:   PetscCall(VecRestoreArray(localF, &farr));
733:   PetscCall(DMRestoreLocalVector(networkdm, &localX));
734:   PetscCall(DMRestoreLocalVector(networkdm, &localXdot));

736:   PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
737:   PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
738:   PetscCall(DMRestoreLocalVector(networkdm, &localF));
739:   PetscFunctionReturn(PETSC_SUCCESS);
740: }

742: /* This function is used for solving the algebraic system only during fault on and
743:    off times. It computes the entire F and then zeros out the part corresponding to
744:    differential equations
745:  F = [0;g(y)];
746: */
747: PetscErrorCode AlgFunction(SNES snes, Vec X, Vec F, void *ctx)
748: {
749:   DM                 networkdm;
750:   Vec                localX, localF;
751:   PetscInt           vfrom, vto, offsetfrom, offsetto;
752:   PetscInt           v, vStart, vEnd, e;
753:   PetscScalar       *farr;
754:   Userctx           *user = (Userctx *)ctx;
755:   const PetscScalar *xarr;
756:   void              *component;
757:   PetscScalar        Vr = 0, Vi = 0;

759:   PetscFunctionBegin;
760:   PetscCall(VecSet(F, 0.0));
761:   PetscCall(SNESGetDM(snes, &networkdm));
762:   PetscCall(DMGetLocalVector(networkdm, &localF));
763:   PetscCall(DMGetLocalVector(networkdm, &localX));
764:   PetscCall(VecSet(localF, 0.0));

766:   /* update ghost values of locaX and locaXdot */
767:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
768:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

770:   PetscCall(VecGetArrayRead(localX, &xarr));
771:   PetscCall(VecGetArray(localF, &farr));

773:   PetscCall(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:     PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghostvtex));
784:     PetscCall(DMNetworkGetNumComponents(networkdm, v, &numComps));

786:     for (j = 0; j < numComps; j++) {
787:       PetscCall(DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL));
788:       if (key == 1) {
789:         PetscInt        nconnedges;
790:         const PetscInt *connedges;

792:         bus = (Bus *)(component);
793:         PetscCall(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:         PetscCall(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:           PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));

821:           Yfti = branch->yft[0];
822:           Yftr = branch->yft[1];

824:           PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
825:           vfrom = cone[0];
826:           vto   = cone[1];

828:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom));
829:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto));

831:           /*From bus and to bus real and imaginary voltages */
832:           Vfr = xarr[offsetfrom];
833:           Vfi = xarr[offsetfrom + 1];
834:           Vtr = xarr[offsetto];
835:           Vti = xarr[offsetto + 1];

837:           if (vfrom == v) {
838:             farr[offsetfrom] += Yftr * Vti + Yfti * Vtr;
839:             farr[offsetfrom + 1] += Yftr * Vtr - Yfti * Vti;
840:           } else {
841:             farr[offsetto] += Yftr * Vfi + Yfti * Vfr;
842:             farr[offsetto + 1] += Yftr * Vfr - Yfti * Vfi;
843:           }
844:         }
845:       } else if (key == 2) {
846:         if (!ghostvtex) {
847:           PetscScalar Eqp, Edp, delta; /* Generator variables */
848:           PetscScalar Id, Iq;          /* Generator dq axis currents */
849:           PetscScalar Vd, Vq, IGr, IGi, Zdq_inv[4], det;
850:           PetscScalar Xdp, Xqp, Rs; /* Generator parameters */

852:           gen = (Gen *)(component);
853:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen));

855:           /* Generator state variables */
856:           Eqp   = xarr[offsetgen];
857:           Edp   = xarr[offsetgen + 1];
858:           delta = xarr[offsetgen + 2];
859:           /* w     = xarr[idx+3]; not being used */
860:           Id = xarr[offsetgen + 4];
861:           Iq = xarr[offsetgen + 5];

863:           /* Generator parameters */
864:           Xdp = gen->Xdp;
865:           Xqp = gen->Xqp;
866:           Rs  = gen->Rs;

868:           /* Set generator differential equation residual functions to zero */
869:           farr[offsetgen]     = 0;
870:           farr[offsetgen + 1] = 0;
871:           farr[offsetgen + 2] = 0;
872:           farr[offsetgen + 3] = 0;

874:           PetscCall(ri2dq(Vr, Vi, delta, &Vd, &Vq));

876:           /* Algebraic equations for stator currents */
877:           det = Rs * Rs + Xdp * Xqp;

879:           Zdq_inv[0] = Rs / det;
880:           Zdq_inv[1] = Xqp / det;
881:           Zdq_inv[2] = -Xdp / det;
882:           Zdq_inv[3] = Rs / det;

884:           farr[offsetgen + 4] = Zdq_inv[0] * (-Edp + Vd) + Zdq_inv[1] * (-Eqp + Vq) + Id;
885:           farr[offsetgen + 5] = Zdq_inv[2] * (-Edp + Vd) + Zdq_inv[3] * (-Eqp + Vq) + Iq;

887:           /* Add generator current injection to network */
888:           PetscCall(dq2ri(Id, Iq, delta, &IGr, &IGi));

890:           farr[offsetbus] -= IGi;
891:           farr[offsetbus + 1] -= IGr;

893:           /* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */
894:         }
895:       } else if (key == 3) {
896:         if (!ghostvtex) {
897:           PetscInt offsetexc;
898:           PetscCall(DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc));
899:           /* Set exciter differential equation residual functions equal to zero*/
900:           farr[offsetexc]     = 0;
901:           farr[offsetexc + 1] = 0;
902:           farr[offsetexc + 2] = 0;
903:         }
904:       } else if (key == 4) {
905:         if (!ghostvtex) {
906:           PetscInt     k, ld_nsegsp, ld_nsegsq;
907:           PetscScalar *ld_alphap, *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;

909:           load = (Load *)(component);

911:           /* Load Parameters */
912:           ld_nsegsp = load->ld_nsegsp;
913:           ld_alphap = load->ld_alphap;
914:           ld_betap  = load->ld_betap;
915:           ld_nsegsq = load->ld_nsegsq;
916:           ld_alphaq = load->ld_alphaq;
917:           ld_betaq  = load->ld_betaq;

919:           PD0 = load->PD0;
920:           QD0 = load->QD0;

922:           Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
923:           Vm2 = Vm * Vm;
924:           Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
925:           PD = QD = 0.0;
926:           for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar((Vm / Vm0), ld_betap[k]);
927:           for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar((Vm / Vm0), ld_betaq[k]);

929:           /* Load currents */
930:           IDr = (PD * Vr + QD * Vi) / Vm2;
931:           IDi = (-QD * Vr + PD * Vi) / Vm2;

933:           farr[offsetbus] += IDi;
934:           farr[offsetbus + 1] += IDr;
935:         }
936:       }
937:     }
938:   }

940:   PetscCall(VecRestoreArrayRead(localX, &xarr));
941:   PetscCall(VecRestoreArray(localF, &farr));
942:   PetscCall(DMRestoreLocalVector(networkdm, &localX));

944:   PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
945:   PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
946:   PetscCall(DMRestoreLocalVector(networkdm, &localF));
947:   PetscFunctionReturn(PETSC_SUCCESS);
948: }

950: int main(int argc, char **argv)
951: {
952:   PetscInt    i, j, *edgelist = NULL, eStart, eEnd, vStart, vEnd;
953:   PetscInt    genj, excj, loadj, componentkey[5];
954:   PetscInt    nc = 1; /* No. of copies (default = 1) */
955:   PetscMPIInt size, rank;
956:   Vec         X, F_alg;
957:   TS          ts;
958:   SNES        snes_alg, snes;
959:   Bus        *bus;
960:   Branch     *branch;
961:   Gen        *gen;
962:   Exc        *exc;
963:   Load       *load;
964:   DM          networkdm;
965: #if defined(PETSC_USE_LOG)
966:   PetscLogStage stage1;
967: #endif
968:   Userctx  user;
969:   KSP      ksp;
970:   PC       pc;
971:   PetscInt numEdges = 0;

973:   PetscFunctionBeginUser;
974:   PetscCall(PetscInitialize(&argc, &argv, "ex9busnetworkops", help));
975:   PetscCall(PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL));
976:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
977:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));

979:   /* Read initial voltage vector and Ybus */
980:   if (rank == 0) PetscCall(read_data(nc, &gen, &exc, &load, &bus, &branch, &edgelist));

982:   PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
983:   PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(Branch), &componentkey[0]));
984:   PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(Bus), &componentkey[1]));
985:   PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(Gen), &componentkey[2]));
986:   PetscCall(DMNetworkRegisterComponent(networkdm, "excstruct", sizeof(Exc), &componentkey[3]));
987:   PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(Load), &componentkey[4]));

989:   PetscCall(PetscLogStageRegister("Create network", &stage1));
990:   PetscCall(PetscLogStagePush(stage1));

992:   /* Set local number of edges and edge connectivity */
993:   if (rank == 0) numEdges = NBRANCH * nc + (nc - 1);
994:   PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
995:   PetscCall(DMNetworkAddSubnetwork(networkdm, NULL, numEdges, edgelist, NULL));

997:   /* Set up the network layout */
998:   PetscCall(DMNetworkLayoutSetUp(networkdm));

1000:   if (rank == 0) PetscCall(PetscFree(edgelist));

1002:   /* Add network components (physical parameters of nodes and branches) and number of variables */
1003:   if (rank == 0) {
1004:     PetscCall(DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd));
1005:     genj  = 0;
1006:     loadj = 0;
1007:     excj  = 0;
1008:     for (i = eStart; i < eEnd; i++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[0], &branch[i - eStart], 0));

1010:     PetscCall(DMNetworkGetVertexRange(networkdm, &vStart, &vEnd));

1012:     for (i = vStart; i < vEnd; i++) {
1013:       PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[1], &bus[i - vStart], 2));
1014:       if (bus[i - vStart].nofgen) {
1015:         for (j = 0; j < bus[i - vStart].nofgen; j++) {
1016:           /* Add generator */
1017:           PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[2], &gen[genj++], 6));
1018:           /* Add exciter */
1019:           PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[3], &exc[excj++], 3));
1020:         }
1021:       }
1022:       if (bus[i - vStart].nofload) {
1023:         for (j = 0; j < bus[i - vStart].nofload; j++) PetscCall(DMNetworkAddComponent(networkdm, i, componentkey[4], &load[loadj++], 0));
1024:       }
1025:     }
1026:   }

1028:   PetscCall(DMSetUp(networkdm));

1030:   if (rank == 0) PetscCall(PetscFree5(bus, gen, load, branch, exc));

1032:   /* for parallel options: Network partitioning and distribution of data */
1033:   if (size > 1) PetscCall(DMNetworkDistribute(&networkdm, 0));
1034:   PetscCall(PetscLogStagePop());

1036:   PetscCall(DMCreateGlobalVector(networkdm, &X));

1038:   PetscCall(SetInitialGuess(networkdm, X));

1040:   /* Options for fault simulation */
1041:   PetscOptionsBegin(PETSC_COMM_WORLD, NULL, "Transient stability fault options", "");
1042:   user.tfaulton  = 0.02;
1043:   user.tfaultoff = 0.05;
1044:   user.Rfault    = 0.0001;
1045:   user.faultbus  = 8;
1046:   PetscCall(PetscOptionsReal("-tfaulton", "", "", user.tfaulton, &user.tfaulton, NULL));
1047:   PetscCall(PetscOptionsReal("-tfaultoff", "", "", user.tfaultoff, &user.tfaultoff, NULL));
1048:   PetscCall(PetscOptionsInt("-faultbus", "", "", user.faultbus, &user.faultbus, NULL));
1049:   user.t0   = 0.0;
1050:   user.tmax = 0.1;
1051:   PetscCall(PetscOptionsReal("-t0", "", "", user.t0, &user.t0, NULL));
1052:   PetscCall(PetscOptionsReal("-tmax", "", "", user.tmax, &user.tmax, NULL));

1054:   PetscCall(PetscMalloc1(18 * nc, &user.ybusfault));
1055:   for (i = 0; i < 18 * nc; i++) user.ybusfault[i] = 0;
1056:   user.ybusfault[user.faultbus * 2 + 1] = 1 / user.Rfault;
1057:   PetscOptionsEnd();

1059:   /* Setup TS solver                                           */
1060:   /*--------------------------------------------------------*/
1061:   PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
1062:   PetscCall(TSSetDM(ts, (DM)networkdm));
1063:   PetscCall(TSSetType(ts, TSCN));

1065:   PetscCall(TSGetSNES(ts, &snes));
1066:   PetscCall(SNESGetKSP(snes, &ksp));
1067:   PetscCall(KSPGetPC(ksp, &pc));
1068:   PetscCall(PCSetType(pc, PCBJACOBI));

1070:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)FormIFunction, &user));
1071:   PetscCall(TSSetMaxTime(ts, user.tfaulton));
1072:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1073:   PetscCall(TSSetTimeStep(ts, 0.01));
1074:   PetscCall(TSSetFromOptions(ts));

1076:   /*user.alg_flg = PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix.
1077:     eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */
1078:   user.alg_flg = PETSC_FALSE;

1080:   /* Prefault period */
1081:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "... (1) Prefault period ... \n"));

1083:   PetscCall(TSSetSolution(ts, X));
1084:   PetscCall(TSSetUp(ts));
1085:   PetscCall(TSSolve(ts, X));

1087:   /* Create the nonlinear solver for solving the algebraic system */
1088:   PetscCall(VecDuplicate(X, &F_alg));

1090:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_alg));
1091:   PetscCall(SNESSetDM(snes_alg, (DM)networkdm));
1092:   PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));
1093:   PetscCall(SNESSetOptionsPrefix(snes_alg, "alg_"));
1094:   PetscCall(SNESSetFromOptions(snes_alg));

1096:   /* Apply disturbance - resistive fault at user.faultbus */
1097:   /* This is done by adding shunt conductance to the diagonal location
1098:      in the Ybus matrix */
1099:   user.alg_flg = PETSC_TRUE;

1101:   /* Solve the algebraic equations */
1102:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (2) Apply disturbance, solve algebraic equations ... \n"));
1103:   PetscCall(SNESSolve(snes_alg, NULL, X));

1105:   /* Disturbance period */
1106:   PetscCall(TSSetTime(ts, user.tfaulton));
1107:   PetscCall(TSSetMaxTime(ts, user.tfaultoff));
1108:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1109:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)FormIFunction, &user));

1111:   user.alg_flg = PETSC_TRUE;
1112:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (3) Disturbance period ... \n"));
1113:   PetscCall(TSSolve(ts, X));

1115:   /* Remove the fault */
1116:   PetscCall(SNESSetFunction(snes_alg, F_alg, AlgFunction, &user));

1118:   user.alg_flg = PETSC_FALSE;
1119:   /* Solve the algebraic equations */
1120:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (4) Remove fault, solve algebraic equations ... \n"));
1121:   PetscCall(SNESSolve(snes_alg, NULL, X));
1122:   PetscCall(SNESDestroy(&snes_alg));

1124:   /* Post-disturbance period */
1125:   PetscCall(TSSetTime(ts, user.tfaultoff));
1126:   PetscCall(TSSetMaxTime(ts, user.tmax));
1127:   PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER));
1128:   PetscCall(TSSetIFunction(ts, NULL, (TSIFunction)FormIFunction, &user));

1130:   user.alg_flg = PETSC_FALSE;
1131:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "\n... (5) Post-disturbance period ... \n"));
1132:   PetscCall(TSSolve(ts, X));

1134:   PetscCall(PetscFree(user.ybusfault));
1135:   PetscCall(VecDestroy(&F_alg));
1136:   PetscCall(VecDestroy(&X));
1137:   PetscCall(DMDestroy(&networkdm));
1138:   PetscCall(TSDestroy(&ts));
1139:   PetscCall(PetscFinalize());
1140:   return 0;
1141: }

1143: /*TEST

1145:    build:
1146:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

1148:    test:
1149:       args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1150:       localrunfiles: X.bin Ybus.bin ex9busnetworkops

1152:    test:
1153:       suffix: 2
1154:       nsize: 2
1155:       args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1156:       localrunfiles: X.bin Ybus.bin ex9busnetworkops

1158: TEST*/