Actual source code: ex9busdmnetwork.c


  2: static char help[] = "This example uses the same problem set up of ex9busdmnetwork.c. \n\
  3: It demonstrates setting and accessing of variables for individual components, instead of \n\
  4: the network vertices (as used in ex9busdmnetwork.c). This is especially useful where vertices \n\
  5: /edges have multiple-components associated with them and one or more components has physics \n\
  6: associated with it. \n\
  7: Input parameters include:\n\
  8:   -nc : number of copies of the base case\n\n";

 10: /*
 11:    This example was modified from ex9busdmnetwork.c.
 12: */

 14: #include <petscts.h>
 15: #include <petscdmnetwork.h>

 17: #define FREQ    60
 18: #define W_S     (2 * PETSC_PI * FREQ)
 19: #define NGEN    3 /* No. of generators in the 9 bus system */
 20: #define NLOAD   3 /* No. of loads in the 9 bus system */
 21: #define NBUS    9 /* No. of buses in the 9 bus system */
 22: #define NBRANCH 9 /* No. of branches in the 9 bus system */

 24: typedef struct {
 25:   PetscInt    id;    /* Bus Number or extended bus name*/
 26:   PetscScalar mbase; /* MVA base of the machine */
 27:   PetscScalar PG;    /* Generator active power output */
 28:   PetscScalar QG;    /* Generator reactive power output */

 30:   /* Generator constants */
 31:   PetscScalar H;    /* Inertia constant */
 32:   PetscScalar Rs;   /* Stator Resistance */
 33:   PetscScalar Xd;   /* d-axis reactance */
 34:   PetscScalar Xdp;  /* d-axis transient reactance */
 35:   PetscScalar Xq;   /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
 36:   PetscScalar Xqp;  /* q-axis transient reactance */
 37:   PetscScalar Td0p; /* d-axis open circuit time constant */
 38:   PetscScalar Tq0p; /* q-axis open circuit time constant */
 39:   PetscScalar M;    /* M = 2*H/W_S */
 40:   PetscScalar D;    /* D = 0.1*M */
 41:   PetscScalar TM;   /* Mechanical Torque */
 42: } Gen;

 44: typedef struct {
 45:   /* Exciter system constants */
 46:   PetscScalar KA;     /* Voltage regulator gain constant */
 47:   PetscScalar TA;     /* Voltage regulator time constant */
 48:   PetscScalar KE;     /* Exciter gain constant */
 49:   PetscScalar TE;     /* Exciter time constant */
 50:   PetscScalar KF;     /* Feedback stabilizer gain constant */
 51:   PetscScalar TF;     /* Feedback stabilizer time constant */
 52:   PetscScalar k1, k2; /* calculating the saturation function SE = k1*exp(k2*Efd) */
 53:   PetscScalar Vref;   /* Voltage regulator voltage setpoint */
 54: } Exc;

 56: typedef struct {
 57:   PetscInt    id;      /* node id */
 58:   PetscInt    nofgen;  /* Number of generators at the bus*/
 59:   PetscInt    nofload; /*  Number of load at the bus*/
 60:   PetscScalar yff[2];  /* yff[0]= imaginary part of admittance, yff[1]=real part of admittance*/
 61:   PetscScalar vr;      /* Real component of bus voltage */
 62:   PetscScalar vi;      /* Imaginary component of bus voltage */
 63: } Bus;

 65: /* Load constants
 66:   We use a composite load model that describes the load and reactive powers at each time instant as follows
 67:   P(t) = \sum\limits_{i=0}^ld_nsegsp \ld_alphap_i*P_D0(\frac{V_m(t)}{V_m0})^\ld_betap_i
 68:   Q(t) = \sum\limits_{i=0}^ld_nsegsq \ld_alphaq_i*Q_D0(\frac{V_m(t)}{V_m0})^\ld_betaq_i
 69:   where
 70:     id                  - index of the load
 71:     ld_nsegsp,ld_nsegsq - Number of individual load models for real and reactive power loads
 72:     ld_alphap,ld_alphap - Percentage contribution (weights) or loads
 73:     P_D0                - Real power load
 74:     Q_D0                - Reactive power load
 75:     Vm(t)              - Voltage magnitude at time t
 76:     Vm0                - Voltage magnitude at t = 0
 77:     ld_betap, ld_betaq  - exponents describing the load model for real and reactive part

 79:     Note: All loads have the same characteristic currently.
 80:   */
 81: typedef struct {
 82:   PetscInt    id; /* bus id */
 83:   PetscInt    ld_nsegsp, ld_nsegsq;
 84:   PetscScalar PD0, QD0;
 85:   PetscScalar ld_alphap[3]; /* ld_alphap=[1,0,0], an array, not a value, so use *ld_alphap; */
 86:   PetscScalar ld_betap[3], ld_alphaq[3], ld_betaq[3];
 87: } Load;

 89: typedef struct {
 90:   PetscInt    id;     /* node id */
 91:   PetscScalar yft[2]; /* yft[0]= imaginary part of admittance, yft[1]=real part of admittance*/
 92: } Branch;

 94: typedef struct {
 95:   PetscReal    tfaulton, tfaultoff; /* Fault on and off times */
 96:   PetscReal    t;
 97:   PetscReal    t0, tmax; /* initial time and final time */
 98:   PetscInt     faultbus; /* Fault bus */
 99:   PetscScalar  Rfault;   /* Fault resistance (pu) */
100:   PetscScalar *ybusfault;
101:   PetscBool    alg_flg;
102: } Userctx;

104: /* Used to read data into the DMNetwork components */
105: PetscErrorCode read_data(PetscInt nc, Gen **pgen, Exc **pexc, Load **pload, Bus **pbus, Branch **pbranch, PetscInt **pedgelist)
106: {
107:   PetscInt           i, j, row[1], col[2];
108:   PetscInt          *edgelist;
109:   PetscInt           nofgen[9]  = {1, 1, 1, 0, 0, 0, 0, 0, 0}; /* Buses at which generators are incident */
110:   PetscInt           nofload[9] = {0, 0, 0, 0, 1, 1, 0, 1, 0}; /* Buses at which loads are incident */
111:   const PetscScalar *varr;
112:   PetscScalar        M[3], D[3];
113:   Bus               *bus;
114:   Branch            *branch;
115:   Gen               *gen;
116:   Exc               *exc;
117:   Load              *load;
118:   Mat                Ybus;
119:   Vec                V0;

121:   /*10 parameters*/
122:   /* Generator real and reactive powers (found via loadflow) */
123:   static const PetscScalar PG[3] = {0.716786142395021, 1.630000000000000, 0.850000000000000};
124:   static const PetscScalar QG[3] = {0.270702180178785, 0.066120127797275, -0.108402221791588};

126:   /* Generator constants */
127:   static const PetscScalar H[3]    = {23.64, 6.4, 3.01};       /* Inertia constant */
128:   static const PetscScalar Rs[3]   = {0.0, 0.0, 0.0};          /* Stator Resistance */
129:   static const PetscScalar Xd[3]   = {0.146, 0.8958, 1.3125};  /* d-axis reactance */
130:   static const PetscScalar Xdp[3]  = {0.0608, 0.1198, 0.1813}; /* d-axis transient reactance */
131:   static const PetscScalar Xq[3]   = {0.4360, 0.8645, 1.2578}; /* q-axis reactance Xq(1) set to 0.4360, value given in text 0.0969 */
132:   static const PetscScalar Xqp[3]  = {0.0969, 0.1969, 0.25};   /* q-axis transient reactance */
133:   static const PetscScalar Td0p[3] = {8.96, 6.0, 5.89};        /* d-axis open circuit time constant */
134:   static const PetscScalar Tq0p[3] = {0.31, 0.535, 0.6};       /* q-axis open circuit time constant */

136:   /* Exciter system constants (8 parameters)*/
137:   static const PetscScalar KA[3] = {20.0, 20.0, 20.0};    /* Voltage regulartor gain constant */
138:   static const PetscScalar TA[3] = {0.2, 0.2, 0.2};       /* Voltage regulator time constant */
139:   static const PetscScalar KE[3] = {1.0, 1.0, 1.0};       /* Exciter gain constant */
140:   static const PetscScalar TE[3] = {0.314, 0.314, 0.314}; /* Exciter time constant */
141:   static const PetscScalar KF[3] = {0.063, 0.063, 0.063}; /* Feedback stabilizer gain constant */
142:   static const PetscScalar TF[3] = {0.35, 0.35, 0.35};    /* Feedback stabilizer time constant */
143:   static const PetscScalar k1[3] = {0.0039, 0.0039, 0.0039};
144:   static const PetscScalar k2[3] = {1.555, 1.555, 1.555}; /* k1 and k2 for calculating the saturation function SE = k1*exp(k2*Efd) */

146:   /* Load constants */
147:   static const PetscScalar PD0[3]       = {1.25, 0.9, 1.0};
148:   static const PetscScalar QD0[3]       = {0.5, 0.3, 0.35};
149:   static const PetscScalar ld_alphaq[3] = {1, 0, 0};
150:   static const PetscScalar ld_betaq[3]  = {2, 1, 0};
151:   static const PetscScalar ld_betap[3]  = {2, 1, 0};
152:   static const PetscScalar ld_alphap[3] = {1, 0, 0};
153:   PetscInt                 ld_nsegsp[3] = {3, 3, 3};
154:   PetscInt                 ld_nsegsq[3] = {3, 3, 3};
155:   PetscViewer              Xview, Ybusview;
156:   PetscInt                 neqs_net, m, n;

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

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

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

173:   MatGetLocalSize(Ybus, &m, &n);
174:   neqs_net = 2 * NBUS; /* # eqs. for network subsystem   */

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

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

187:   VecGetArrayRead(V0, &varr);

189:   /* read bus data */
190:   for (i = 0; i < nc; i++) {
191:     for (j = 0; j < NBUS; j++) {
192:       bus[i * 9 + j].id      = i * 9 + j;
193:       bus[i * 9 + j].nofgen  = nofgen[j];
194:       bus[i * 9 + j].nofload = nofload[j];
195:       bus[i * 9 + j].vr      = varr[2 * j];
196:       bus[i * 9 + j].vi      = varr[2 * j + 1];
197:       row[0]                 = 2 * j;
198:       col[0]                 = 2 * j;
199:       col[1]                 = 2 * j + 1;
200:       /* real and imaginary part of admittance from Ybus into yff */
201:       MatGetValues(Ybus, 1, row, 2, col, bus[i * 9 + j].yff);
202:     }
203:   }

205:   /* read generator data */
206:   for (i = 0; i < nc; i++) {
207:     for (j = 0; j < NGEN; j++) {
208:       gen[i * 3 + j].id   = i * 3 + j;
209:       gen[i * 3 + j].PG   = PG[j];
210:       gen[i * 3 + j].QG   = QG[j];
211:       gen[i * 3 + j].H    = H[j];
212:       gen[i * 3 + j].Rs   = Rs[j];
213:       gen[i * 3 + j].Xd   = Xd[j];
214:       gen[i * 3 + j].Xdp  = Xdp[j];
215:       gen[i * 3 + j].Xq   = Xq[j];
216:       gen[i * 3 + j].Xqp  = Xqp[j];
217:       gen[i * 3 + j].Td0p = Td0p[j];
218:       gen[i * 3 + j].Tq0p = Tq0p[j];
219:       gen[i * 3 + j].M    = M[j];
220:       gen[i * 3 + j].D    = D[j];
221:     }
222:   }

224:   for (i = 0; i < nc; i++) {
225:     for (j = 0; j < NGEN; j++) {
226:       /* exciter system */
227:       exc[i * 3 + j].KA = KA[j];
228:       exc[i * 3 + j].TA = TA[j];
229:       exc[i * 3 + j].KE = KE[j];
230:       exc[i * 3 + j].TE = TE[j];
231:       exc[i * 3 + j].KF = KF[j];
232:       exc[i * 3 + j].TF = TF[j];
233:       exc[i * 3 + j].k1 = k1[j];
234:       exc[i * 3 + j].k2 = k2[j];
235:     }
236:   }

238:   /* read load data */
239:   for (i = 0; i < nc; i++) {
240:     for (j = 0; j < NLOAD; j++) {
241:       load[i * 3 + j].id        = i * 3 + j;
242:       load[i * 3 + j].PD0       = PD0[j];
243:       load[i * 3 + j].QD0       = QD0[j];
244:       load[i * 3 + j].ld_nsegsp = ld_nsegsp[j];

246:       load[i * 3 + j].ld_alphap[0] = ld_alphap[0];
247:       load[i * 3 + j].ld_alphap[1] = ld_alphap[1];
248:       load[i * 3 + j].ld_alphap[2] = ld_alphap[2];

250:       load[i * 3 + j].ld_alphaq[0] = ld_alphaq[0];
251:       load[i * 3 + j].ld_alphaq[1] = ld_alphaq[1];
252:       load[i * 3 + j].ld_alphaq[2] = ld_alphaq[2];

254:       load[i * 3 + j].ld_betap[0] = ld_betap[0];
255:       load[i * 3 + j].ld_betap[1] = ld_betap[1];
256:       load[i * 3 + j].ld_betap[2] = ld_betap[2];
257:       load[i * 3 + j].ld_nsegsq   = ld_nsegsq[j];

259:       load[i * 3 + j].ld_betaq[0] = ld_betaq[0];
260:       load[i * 3 + j].ld_betaq[1] = ld_betaq[1];
261:       load[i * 3 + j].ld_betaq[2] = ld_betaq[2];
262:     }
263:   }
264:   PetscCalloc1(2 * NBRANCH * nc + 2 * (nc - 1), &edgelist);

266:   /* read edgelist */
267:   for (i = 0; i < nc; i++) {
268:     for (j = 0; j < NBRANCH; j++) {
269:       switch (j) {
270:       case 0:
271:         edgelist[i * 18 + 2 * j]     = 0 + 9 * i;
272:         edgelist[i * 18 + 2 * j + 1] = 3 + 9 * i;
273:         break;
274:       case 1:
275:         edgelist[i * 18 + 2 * j]     = 1 + 9 * i;
276:         edgelist[i * 18 + 2 * j + 1] = 6 + 9 * i;
277:         break;
278:       case 2:
279:         edgelist[i * 18 + 2 * j]     = 2 + 9 * i;
280:         edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
281:         break;
282:       case 3:
283:         edgelist[i * 18 + 2 * j]     = 3 + 9 * i;
284:         edgelist[i * 18 + 2 * j + 1] = 4 + 9 * i;
285:         break;
286:       case 4:
287:         edgelist[i * 18 + 2 * j]     = 3 + 9 * i;
288:         edgelist[i * 18 + 2 * j + 1] = 5 + 9 * i;
289:         break;
290:       case 5:
291:         edgelist[i * 18 + 2 * j]     = 4 + 9 * i;
292:         edgelist[i * 18 + 2 * j + 1] = 6 + 9 * i;
293:         break;
294:       case 6:
295:         edgelist[i * 18 + 2 * j]     = 5 + 9 * i;
296:         edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
297:         break;
298:       case 7:
299:         edgelist[i * 18 + 2 * j]     = 6 + 9 * i;
300:         edgelist[i * 18 + 2 * j + 1] = 7 + 9 * i;
301:         break;
302:       case 8:
303:         edgelist[i * 18 + 2 * j]     = 7 + 9 * i;
304:         edgelist[i * 18 + 2 * j + 1] = 8 + 9 * i;
305:         break;
306:       default:
307:         break;
308:       }
309:     }
310:   }

312:   /* for connecting last bus of previous network(9*i-1) to first bus of next network(9*i), the branch admittance=-0.0301407+j17.3611 */
313:   for (i = 1; i < nc; i++) {
314:     edgelist[18 * nc + 2 * (i - 1)]     = 8 + (i - 1) * 9;
315:     edgelist[18 * nc + 2 * (i - 1) + 1] = 9 * i;

317:     /* adding admittances to the off-diagonal elements */
318:     branch[9 * nc + (i - 1)].id     = 9 * nc + (i - 1);
319:     branch[9 * nc + (i - 1)].yft[0] = 17.3611;
320:     branch[9 * nc + (i - 1)].yft[1] = -0.0301407;

322:     /* subtracting admittances from the diagonal elements */
323:     bus[9 * i - 1].yff[0] -= 17.3611;
324:     bus[9 * i - 1].yff[1] -= -0.0301407;

326:     bus[9 * i].yff[0] -= 17.3611;
327:     bus[9 * i].yff[1] -= -0.0301407;
328:   }

330:   /* read branch data */
331:   for (i = 0; i < nc; i++) {
332:     for (j = 0; j < NBRANCH; j++) {
333:       branch[i * 9 + j].id = i * 9 + j;

335:       row[0] = edgelist[2 * j] * 2;
336:       col[0] = edgelist[2 * j + 1] * 2;
337:       col[1] = edgelist[2 * j + 1] * 2 + 1;
338:       MatGetValues(Ybus, 1, row, 2, col, branch[i * 9 + j].yft); /*imaginary part of admittance*/
339:     }
340:   }

342:   *pgen      = gen;
343:   *pexc      = exc;
344:   *pload     = load;
345:   *pbus      = bus;
346:   *pbranch   = branch;
347:   *pedgelist = edgelist;

349:   VecRestoreArrayRead(V0, &varr);

351:   /* Destroy unnecessary stuff */
352:   MatDestroy(&Ybus);
353:   VecDestroy(&V0);
354:   return 0;
355: }

357: PetscErrorCode SetInitialGuess(DM networkdm, Vec X)
358: {
359:   Bus         *bus;
360:   Gen         *gen;
361:   Exc         *exc;
362:   PetscInt     v, vStart, vEnd, offset;
363:   PetscInt     key, numComps, j;
364:   PetscBool    ghostvtex;
365:   Vec          localX;
366:   PetscScalar *xarr;
367:   PetscScalar  Vr = 0, Vi = 0, Vm = 0, Vm2; /* Terminal voltage variables */
368:   PetscScalar  IGr, IGi;                    /* Generator real and imaginary current */
369:   PetscScalar  Eqp, Edp, delta;             /* Generator variables */
370:   PetscScalar  Efd = 0, RF, VR;             /* Exciter variables */
371:   PetscScalar  Vd, Vq;                      /* Generator dq axis voltages */
372:   PetscScalar  Id, Iq;                      /* Generator dq axis currents */
373:   PetscScalar  theta;                       /* Generator phase angle */
374:   PetscScalar  SE;
375:   void        *component;

377:   DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);
378:   DMGetLocalVector(networkdm, &localX);

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

384:   VecGetArray(localX, &xarr);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

473: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
474: PetscErrorCode ri2dq(PetscScalar Fr, PetscScalar Fi, PetscScalar delta, PetscScalar *Fd, PetscScalar *Fq)
475: {
476:   *Fd = Fr * PetscSinScalar(delta) - Fi * PetscCosScalar(delta);
477:   *Fq = Fr * PetscCosScalar(delta) + Fi * PetscSinScalar(delta);
478:   return 0;
479: }

481: /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
482: PetscErrorCode FormIFunction(TS ts, PetscReal t, Vec X, Vec Xdot, Vec F, Userctx *user)
483: {
484:   DM                 networkdm;
485:   Vec                localX, localXdot, localF;
486:   PetscInt           vfrom, vto, offsetfrom, offsetto;
487:   PetscInt           v, vStart, vEnd, e;
488:   PetscScalar       *farr;
489:   PetscScalar        Vd = 0, Vq = 0, SE;
490:   const PetscScalar *xarr, *xdotarr;
491:   void              *component;
492:   PetscScalar        Vr = 0, Vi = 0;

494:   VecSet(F, 0.0);

496:   TSGetDM(ts, &networkdm);
497:   DMGetLocalVector(networkdm, &localF);
498:   DMGetLocalVector(networkdm, &localX);
499:   DMGetLocalVector(networkdm, &localXdot);
500:   VecSet(localF, 0.0);

502:   /* update ghost values of localX and localXdot */
503:   DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
504:   DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);

506:   DMGlobalToLocalBegin(networkdm, Xdot, INSERT_VALUES, localXdot);
507:   DMGlobalToLocalEnd(networkdm, Xdot, INSERT_VALUES, localXdot);

509:   VecGetArrayRead(localX, &xarr);
510:   VecGetArrayRead(localXdot, &xdotarr);
511:   VecGetArray(localF, &farr);

513:   DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);

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

528:     DMNetworkIsGhostVertex(networkdm, v, &ghostvtex);
529:     DMNetworkGetNumComponents(networkdm, v, &numComps);

531:     for (j = 0; j < numComps; j++) {
532:       DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL);
533:       if (key == 1) {
534:         PetscInt        nconnedges;
535:         const PetscInt *connedges;

537:         bus = (Bus *)(component);
538:         DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus);
539:         if (!ghostvtex) {
540:           Vr = xarr[offsetbus];
541:           Vi = xarr[offsetbus + 1];

543:           Yffr = bus->yff[1];
544:           Yffi = bus->yff[0];

546:           if (user->alg_flg) {
547:             Yffr += user->ybusfault[bus->id * 2 + 1];
548:             Yffi += user->ybusfault[bus->id * 2];
549:           }
550:           Vr0 = bus->vr;
551:           Vi0 = bus->vi;

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

560:         DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges);

562:         for (i = 0; i < nconnedges; i++) {
563:           Branch         *branch;
564:           PetscInt        keye;
565:           PetscScalar     Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
566:           const PetscInt *cone;

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

571:           Yfti = branch->yft[0];
572:           Yftr = branch->yft[1];

574:           DMNetworkGetConnectedVertices(networkdm, e, &cone);

576:           vfrom = cone[0];
577:           vto   = cone[1];

579:           DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom);
580:           DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto);

582:           /* From bus and to bus real and imaginary voltages */
583:           Vfr = xarr[offsetfrom];
584:           Vfi = xarr[offsetfrom + 1];
585:           Vtr = xarr[offsetto];
586:           Vti = xarr[offsetto + 1];

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

604:           gen = (Gen *)(component);
605:           DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen);

607:           /* Generator state variables */
608:           Eqp   = xarr[offsetgen];
609:           Edp   = xarr[offsetgen + 1];
610:           delta = xarr[offsetgen + 2];
611:           w     = xarr[offsetgen + 3];
612:           Id    = xarr[offsetgen + 4];
613:           Iq    = xarr[offsetgen + 5];

615:           /* Generator parameters */
616:           Xd   = gen->Xd;
617:           Xdp  = gen->Xdp;
618:           Td0p = gen->Td0p;
619:           Xq   = gen->Xq;
620:           Xqp  = gen->Xqp;
621:           Tq0p = gen->Tq0p;
622:           TM   = gen->TM;
623:           D    = gen->D;
624:           M    = gen->M;
625:           Rs   = gen->Rs;

627:           DMNetworkGetLocalVecOffset(networkdm, v, 2, &offsetexc);
628:           Efd = xarr[offsetexc];

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

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

638:           /* Algebraic equations for stator currents */
639:           det = Rs * Rs + Xdp * Xqp;

641:           Zdq_inv[0] = Rs / det;
642:           Zdq_inv[1] = Xqp / det;
643:           Zdq_inv[2] = -Xdp / det;
644:           Zdq_inv[3] = Rs / det;

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

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

651:           /* Add generator current injection to network */
652:           farr[offsetbus] -= IGi;
653:           farr[offsetbus + 1] -= IGr;
654:         }
655:       } else if (key == 3) {
656:         if (!ghostvtex) {
657:           PetscScalar k1, k2, KE, TE, TF, KA, KF, Vref, TA; /* Generator parameters */
658:           PetscScalar Efd, RF, VR;                          /* Exciter variables */

660:           exc = (Exc *)(component);
661:           DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetexc);

663:           Efd = xarr[offsetexc];
664:           RF  = xarr[offsetexc + 1];
665:           VR  = xarr[offsetexc + 2];

667:           k1   = exc->k1;
668:           k2   = exc->k2;
669:           KE   = exc->KE;
670:           TE   = exc->TE;
671:           TF   = exc->TF;
672:           KA   = exc->KA;
673:           KF   = exc->KF;
674:           Vref = exc->Vref;
675:           TA   = exc->TA;

677:           Vm = PetscSqrtScalar(Vd * Vd + Vq * Vq);
678:           SE = k1 * PetscExpScalar(k2 * Efd);

680:           /* Exciter differential equations */
681:           farr[offsetexc]     = (KE * Efd + SE - VR) / TE + xdotarr[offsetexc];
682:           farr[offsetexc + 1] = (RF - KF * Efd / TF) / TF + xdotarr[offsetexc + 1];
683:           farr[offsetexc + 2] = (VR - KA * RF + KA * KF * Efd / TF - KA * (Vref - Vm)) / TA + xdotarr[offsetexc + 2];
684:         }
685:       } else if (key == 4) {
686:         if (!ghostvtex) {
687:           PetscInt     k;
688:           PetscInt     ld_nsegsp;
689:           PetscInt     ld_nsegsq;
690:           PetscScalar *ld_alphap;
691:           PetscScalar *ld_betap, *ld_alphaq, *ld_betaq, PD0, QD0, IDr, IDi;

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

695:           /* Load Parameters */
696:           ld_nsegsp = load->ld_nsegsp;
697:           ld_alphap = load->ld_alphap;
698:           ld_betap  = load->ld_betap;
699:           ld_nsegsq = load->ld_nsegsq;
700:           ld_alphaq = load->ld_alphaq;
701:           ld_betaq  = load->ld_betaq;
702:           PD0       = load->PD0;
703:           QD0       = load->QD0;

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

714:           /* Load currents */
715:           IDr = (PD * Vr + QD * Vi) / Vm2;
716:           IDi = (-QD * Vr + PD * Vi) / Vm2;

718:           /* Load current contribution to the network */
719:           farr[offsetbus] += IDi;
720:           farr[offsetbus + 1] += IDr;
721:         }
722:       }
723:     }
724:   }

726:   VecRestoreArrayRead(localX, &xarr);
727:   VecRestoreArrayRead(localXdot, &xdotarr);
728:   VecRestoreArray(localF, &farr);
729:   DMRestoreLocalVector(networkdm, &localX);
730:   DMRestoreLocalVector(networkdm, &localXdot);

732:   DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F);
733:   DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F);
734:   DMRestoreLocalVector(networkdm, &localF);
735:   return 0;
736: }

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

755:   VecSet(F, 0.0);
756:   SNESGetDM(snes, &networkdm);
757:   DMGetLocalVector(networkdm, &localF);
758:   DMGetLocalVector(networkdm, &localX);
759:   VecSet(localF, 0.0);

761:   /* update ghost values of locaX and locaXdot */
762:   DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
763:   DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);

765:   VecGetArrayRead(localX, &xarr);
766:   VecGetArray(localF, &farr);

768:   DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);

770:   for (v = vStart; v < vEnd; v++) {
771:     PetscInt    i, j, offsetbus, offsetgen, key, numComps;
772:     PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0 = 0, Vi0 = 0, PD, QD;
773:     Bus        *bus;
774:     Gen        *gen;
775:     Load       *load;
776:     PetscBool   ghostvtex;

778:     DMNetworkIsGhostVertex(networkdm, v, &ghostvtex);
779:     DMNetworkGetNumComponents(networkdm, v, &numComps);

781:     for (j = 0; j < numComps; j++) {
782:       DMNetworkGetComponent(networkdm, v, j, &key, &component, NULL);
783:       if (key == 1) {
784:         PetscInt        nconnedges;
785:         const PetscInt *connedges;

787:         bus = (Bus *)(component);
788:         DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetbus);
789:         if (!ghostvtex) {
790:           Vr = xarr[offsetbus];
791:           Vi = xarr[offsetbus + 1];

793:           Yffr = bus->yff[1];
794:           Yffi = bus->yff[0];
795:           if (user->alg_flg) {
796:             Yffr += user->ybusfault[bus->id * 2 + 1];
797:             Yffi += user->ybusfault[bus->id * 2];
798:           }
799:           Vr0 = bus->vr;
800:           Vi0 = bus->vi;

802:           farr[offsetbus] += Yffi * Vr + Yffr * Vi;
803:           farr[offsetbus + 1] += Yffr * Vr - Yffi * Vi;
804:         }
805:         DMNetworkGetSupportingEdges(networkdm, v, &nconnedges, &connedges);

807:         for (i = 0; i < nconnedges; i++) {
808:           Branch         *branch;
809:           PetscInt        keye;
810:           PetscScalar     Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
811:           const PetscInt *cone;

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

816:           Yfti = branch->yft[0];
817:           Yftr = branch->yft[1];

819:           DMNetworkGetConnectedVertices(networkdm, e, &cone);
820:           vfrom = cone[0];
821:           vto   = cone[1];

823:           DMNetworkGetLocalVecOffset(networkdm, vfrom, 0, &offsetfrom);
824:           DMNetworkGetLocalVecOffset(networkdm, vto, 0, &offsetto);

826:           /*From bus and to bus real and imaginary voltages */
827:           Vfr = xarr[offsetfrom];
828:           Vfi = xarr[offsetfrom + 1];
829:           Vtr = xarr[offsetto];
830:           Vti = xarr[offsetto + 1];

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

847:           gen = (Gen *)(component);
848:           DMNetworkGetLocalVecOffset(networkdm, v, j, &offsetgen);

850:           /* Generator state variables */
851:           Eqp   = xarr[offsetgen];
852:           Edp   = xarr[offsetgen + 1];
853:           delta = xarr[offsetgen + 2];
854:           /* w     = xarr[idx+3]; not being used */
855:           Id = xarr[offsetgen + 4];
856:           Iq = xarr[offsetgen + 5];

858:           /* Generator parameters */
859:           Xdp = gen->Xdp;
860:           Xqp = gen->Xqp;
861:           Rs  = gen->Rs;

863:           /* Set generator differential equation residual functions to zero */
864:           farr[offsetgen]     = 0;
865:           farr[offsetgen + 1] = 0;
866:           farr[offsetgen + 2] = 0;
867:           farr[offsetgen + 3] = 0;

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

871:           /* Algebraic equations for stator currents */
872:           det = Rs * Rs + Xdp * Xqp;

874:           Zdq_inv[0] = Rs / det;
875:           Zdq_inv[1] = Xqp / det;
876:           Zdq_inv[2] = -Xdp / det;
877:           Zdq_inv[3] = Rs / det;

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

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

885:           farr[offsetbus] -= IGi;
886:           farr[offsetbus + 1] -= IGr;

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

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

906:           /* Load Parameters */
907:           ld_nsegsp = load->ld_nsegsp;
908:           ld_alphap = load->ld_alphap;
909:           ld_betap  = load->ld_betap;
910:           ld_nsegsq = load->ld_nsegsq;
911:           ld_alphaq = load->ld_alphaq;
912:           ld_betaq  = load->ld_betaq;

914:           PD0 = load->PD0;
915:           QD0 = load->QD0;

917:           Vm  = PetscSqrtScalar(Vr * Vr + Vi * Vi);
918:           Vm2 = Vm * Vm;
919:           Vm0 = PetscSqrtScalar(Vr0 * Vr0 + Vi0 * Vi0);
920:           PD = QD = 0.0;
921:           for (k = 0; k < ld_nsegsp; k++) PD += ld_alphap[k] * PD0 * PetscPowScalar((Vm / Vm0), ld_betap[k]);
922:           for (k = 0; k < ld_nsegsq; k++) QD += ld_alphaq[k] * QD0 * PetscPowScalar((Vm / Vm0), ld_betaq[k]);

924:           /* Load currents */
925:           IDr = (PD * Vr + QD * Vi) / Vm2;
926:           IDi = (-QD * Vr + PD * Vi) / Vm2;

928:           farr[offsetbus] += IDi;
929:           farr[offsetbus + 1] += IDr;
930:         }
931:       }
932:     }
933:   }

935:   VecRestoreArrayRead(localX, &xarr);
936:   VecRestoreArray(localF, &farr);
937:   DMRestoreLocalVector(networkdm, &localX);

939:   DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F);
940:   DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F);
941:   DMRestoreLocalVector(networkdm, &localF);
942:   return 0;
943: }

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

969:   PetscInitialize(&argc, &argv, "ex9busnetworkops", help);
970:   PetscOptionsGetInt(NULL, NULL, "-nc", &nc, NULL);
971:   MPI_Comm_size(PETSC_COMM_WORLD, &size);
972:   MPI_Comm_rank(PETSC_COMM_WORLD, &rank);

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

977:   DMNetworkCreate(PETSC_COMM_WORLD, &networkdm);
978:   DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(Branch), &componentkey[0]);
979:   DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(Bus), &componentkey[1]);
980:   DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(Gen), &componentkey[2]);
981:   DMNetworkRegisterComponent(networkdm, "excstruct", sizeof(Exc), &componentkey[3]);
982:   DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(Load), &componentkey[4]);

984:   PetscLogStageRegister("Create network", &stage1);
985:   PetscLogStagePush(stage1);

987:   /* Set local number of edges and edge connectivity */
988:   if (rank == 0) numEdges = NBRANCH * nc + (nc - 1);
989:   DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1);
990:   DMNetworkAddSubnetwork(networkdm, NULL, numEdges, edgelist, NULL);

992:   /* Set up the network layout */
993:   DMNetworkLayoutSetUp(networkdm);

995:   if (rank == 0) PetscFree(edgelist);

997:   /* Add network components (physical parameters of nodes and branches) and number of variables */
998:   if (rank == 0) {
999:     DMNetworkGetEdgeRange(networkdm, &eStart, &eEnd);
1000:     genj  = 0;
1001:     loadj = 0;
1002:     excj  = 0;
1003:     for (i = eStart; i < eEnd; i++) DMNetworkAddComponent(networkdm, i, componentkey[0], &branch[i - eStart], 0);

1005:     DMNetworkGetVertexRange(networkdm, &vStart, &vEnd);

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

1023:   DMSetUp(networkdm);

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

1027:   /* for parallel options: Network partitioning and distribution of data */
1028:   if (size > 1) DMNetworkDistribute(&networkdm, 0);
1029:   PetscLogStagePop();

1031:   DMCreateGlobalVector(networkdm, &X);

1033:   SetInitialGuess(networkdm, X);

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

1049:   PetscMalloc1(18 * nc, &user.ybusfault);
1050:   for (i = 0; i < 18 * nc; i++) user.ybusfault[i] = 0;
1051:   user.ybusfault[user.faultbus * 2 + 1] = 1 / user.Rfault;
1052:   PetscOptionsEnd();

1054:   /* Setup TS solver                                           */
1055:   /*--------------------------------------------------------*/
1056:   TSCreate(PETSC_COMM_WORLD, &ts);
1057:   TSSetDM(ts, (DM)networkdm);
1058:   TSSetType(ts, TSCN);

1060:   TSGetSNES(ts, &snes);
1061:   SNESGetKSP(snes, &ksp);
1062:   KSPGetPC(ksp, &pc);
1063:   PCSetType(pc, PCBJACOBI);

1065:   TSSetIFunction(ts, NULL, (TSIFunction)FormIFunction, &user);
1066:   TSSetMaxTime(ts, user.tfaulton);
1067:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1068:   TSSetTimeStep(ts, 0.01);
1069:   TSSetFromOptions(ts);

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

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

1078:   TSSetSolution(ts, X);
1079:   TSSetUp(ts);
1080:   TSSolve(ts, X);

1082:   /* Create the nonlinear solver for solving the algebraic system */
1083:   VecDuplicate(X, &F_alg);

1085:   SNESCreate(PETSC_COMM_WORLD, &snes_alg);
1086:   SNESSetDM(snes_alg, (DM)networkdm);
1087:   SNESSetFunction(snes_alg, F_alg, AlgFunction, &user);
1088:   SNESSetOptionsPrefix(snes_alg, "alg_");
1089:   SNESSetFromOptions(snes_alg);

1091:   /* Apply disturbance - resistive fault at user.faultbus */
1092:   /* This is done by adding shunt conductance to the diagonal location
1093:      in the Ybus matrix */
1094:   user.alg_flg = PETSC_TRUE;

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

1100:   /* Disturbance period */
1101:   TSSetTime(ts, user.tfaulton);
1102:   TSSetMaxTime(ts, user.tfaultoff);
1103:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1104:   TSSetIFunction(ts, NULL, (TSIFunction)FormIFunction, &user);

1106:   user.alg_flg = PETSC_TRUE;
1107:   if (rank == 0) PetscPrintf(PETSC_COMM_SELF, "\n... (3) Disturbance period ... \n");
1108:   TSSolve(ts, X);

1110:   /* Remove the fault */
1111:   SNESSetFunction(snes_alg, F_alg, AlgFunction, &user);

1113:   user.alg_flg = PETSC_FALSE;
1114:   /* Solve the algebraic equations */
1115:   if (rank == 0) PetscPrintf(PETSC_COMM_SELF, "\n... (4) Remove fault, solve algebraic equations ... \n");
1116:   SNESSolve(snes_alg, NULL, X);
1117:   SNESDestroy(&snes_alg);

1119:   /* Post-disturbance period */
1120:   TSSetTime(ts, user.tfaultoff);
1121:   TSSetMaxTime(ts, user.tmax);
1122:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
1123:   TSSetIFunction(ts, NULL, (TSIFunction)FormIFunction, &user);

1125:   user.alg_flg = PETSC_FALSE;
1126:   if (rank == 0) PetscPrintf(PETSC_COMM_SELF, "\n... (5) Post-disturbance period ... \n");
1127:   TSSolve(ts, X);

1129:   PetscFree(user.ybusfault);
1130:   VecDestroy(&F_alg);
1131:   VecDestroy(&X);
1132:   DMDestroy(&networkdm);
1133:   TSDestroy(&ts);
1134:   PetscFinalize();
1135:   return 0;
1136: }

1138: /*TEST

1140:    build:
1141:       requires: double !complex !defined(PETSC_USE_64BIT_INDICES)

1143:    test:
1144:       args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1145:       localrunfiles: X.bin Ybus.bin ex9busnetworkops

1147:    test:
1148:       suffix: 2
1149:       nsize: 2
1150:       args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
1151:       localrunfiles: X.bin Ybus.bin ex9busnetworkops

1153: TEST*/