2: static char help[] = "Power grid stability analysis of a large power system.\n\
  3: The base case is the WECC 9 bus system based on the example given in the book Power\n\
  4: Systems Dynamics and Stability (Chapter 7) by P. Sauer and M. A. Pai.\n\
  5: The base power grid in this example consists of 9 buses (nodes), 3 generators,\n\
  6: 3 loads, and 9 transmission lines. The network equations are written\n\
  7: in current balance form using rectangular coordiantes. \n\
  8: DMNetwork is used to manage the variables and equations of the entire system.\n\
  9: Input parameters include:\n\
 10:   -nc : number of copies of the base case\n\n";

 12: /* T
 13:    Concepts: DMNetwork
 14:    Concepts: PETSc TS solver

 16:    This example was contributed by Bikiran Guha and Jianqiao Huang, Illinois Institute of Technology, 2017.
 17: */

 19: #include <petscts.h>
 20: #include <petscdmnetwork.h>

 22: #define FREQ 60
 23: #define W_S (2*PETSC_PI*FREQ)
 24: #define NGEN    3   /* No. of generators in the 9 bus system */
 25: #define NLOAD   3   /* No. of loads in the 9 bus system */
 26: #define NBUS    9   /* No. of buses in the 9 bus system */
 27: #define NBRANCH 9   /* No. of branches in the 9 bus system */

 29: typedef struct {
 30:   PetscInt    id;    /* Bus Number or extended bus name*/
 31:   PetscScalar mbase; /* MVA base of the machine */
 32:   PetscScalar PG;    /* Generator active power output */
 33:   PetscScalar QG;    /* Generator reactive power output */

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

 48:   /* Exciter system constants */
 49:   PetscScalar KA ;   /* Voltage regulator gain constant */
 50:   PetscScalar TA;    /* Voltage regulator time constant */
 51:   PetscScalar KE;    /* Exciter gain constant */
 52:   PetscScalar TE;    /* Exciter time constant */
 53:   PetscScalar KF;    /* Feedback stabilizer gain constant */
 54:   PetscScalar TF;    /* Feedback stabilizer time constant */
 55:   PetscScalar k1,k2; /* calculating the saturation function SE = k1*exp(k2*Efd) */
 56:   PetscScalar Vref;  /* Voltage regulator voltage setpoint */
 57: } Gen;

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

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

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

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

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

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

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

128:   /* Generator constants */
129:   static const PetscScalar H[3]    = {23.64,6.4,3.01};   /* Inertia constant */
130:   static const PetscScalar Rs[3]   = {0.0,0.0,0.0}; /* Stator Resistance */
131:   static const PetscScalar Xd[3]   = {0.146,0.8958,1.3125};  /* d-axis reactance */
132:   static const PetscScalar Xdp[3]  = {0.0608,0.1198,0.1813}; /* d-axis transient reactance */
133:   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 */
134:   static const PetscScalar Xqp[3]  = {0.0969,0.1969,0.25}; /* q-axis transient reactance */
135:   static const PetscScalar Td0p[3] = {8.96,6.0,5.89}; /* d-axis open circuit time constant */
136:   static const PetscScalar Tq0p[3] = {0.31,0.535,0.6}; /* q-axis open circuit time constant */

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

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

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

167:    MatCreate(PETSC_COMM_SELF,&Ybus);
168:    MatSetType(Ybus,MATBAIJ);
169:    MatLoad(Ybus,Ybusview);

171:    /* Destroy unnecessary stuff */
172:    PetscViewerDestroy(&Xview);
173:    PetscViewerDestroy(&Ybusview);

175:    MatGetLocalSize(Ybus,&m,&n);
176:    neqs_net = 2*NBUS; /* # eqs. for network subsystem   */
177:    if (m != neqs_net || n != neqs_net) SETERRQ(PETSC_COMM_SELF,0,"matrix Ybus is in wrong sizes");

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

186:    /* Alocate memory for total number of buses, generators, loads and branches */
187:    PetscCalloc4(NBUS*nc,&bus,NGEN*nc,&gen,NLOAD*nc,&load,NBRANCH*nc+(nc-1),&branch);

189:    VecGetArray(V0,&varr);

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

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

224:        /* exciter system */
225:        gen[i*3+j].KA = KA[j];
226:        gen[i*3+j].TA = TA[j];
227:        gen[i*3+j].KE = KE[j];
228:        gen[i*3+j].TE = TE[j];
229:        gen[i*3+j].KF = KF[j];
230:        gen[i*3+j].TF = TF[j];
231:        gen[i*3+j].k1 = k1[j];
232:        gen[i*3+j].k2 = k2[j];
233:      }
234:    }

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

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

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

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

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

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

310:    /* 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 */
311:     for (i = 1; i<nc; i++){
312:         edgelist[18*nc+2*(i-1)]   = 8+(i-1)*9;
313:         edgelist[18*nc+2*(i-1)+1] = 9*i;

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

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

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

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

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

340:    *pgen      = gen;
341:    *pload     = load;
342:    *pbus      = bus;
343:    *pbranch   = branch;
344:    *pedgelist = edgelist;

346:    VecRestoreArray(V0,&varr);

348:    /* Destroy unnecessary stuff */
349:    MatDestroy(&Ybus);
350:    VecDestroy(&V0);
351:    return(0);
352: }

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

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

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

383:   VecGetArray(localX,&xarr);
384:   DMNetworkGetComponentDataArray(networkdm,&arr);

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

390:     DMNetworkGetVariableOffset(networkdm,v,&offset);
391:     DMNetworkGetNumComponents(networkdm,v,&numComps);
392:     for (j=0; j < numComps; j++) {
393:       DMNetworkGetComponentKeyOffset(networkdm,v,j,&key,&offsetd);
394:       if (key == 1) {
395:         bus = (Bus*)(arr+offsetd);

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*)(arr+offsetd);

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;
432:         idx     = offset+2;

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

441:         /* Exciter */
442:         Efd = Eqp + (gen->Xd - gen->Xdp)*Id;
443:         SE  = gen->k1*PetscExpScalar(gen->k2*Efd);
444:         VR  = gen->KE*Efd + SE;
445:         RF  = gen->KF*Efd/gen->TF;

447:         xarr[idx+6] = Efd;
448:         xarr[idx+7] = RF;
449:         xarr[idx+8] = VR;

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

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

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

480: /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
481: PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx *user)
482: {
483:   PetscErrorCode                     ierr;
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,Vq,SE;
490:   const PetscScalar                  *xarr,*xdotarr;
491:   DMNetworkComponentGenericDataType  *arr;

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);
514:   DMNetworkGetComponentDataArray(networkdm,&arr);

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

529:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
530:     DMNetworkGetNumComponents(networkdm,v,&numComps);
531:     DMNetworkGetVariableOffset(networkdm,v,&offset);

533:     for (j = 0; j < numComps; j++) {
534:       DMNetworkGetComponentKeyOffset(networkdm,v,j,&key,&offsetd);
535:       if (key == 1) {
536:         PetscInt       nconnedges;
537:         const PetscInt *connedges;

539:         bus = (Bus*)(arr+offsetd);
540:         if (!ghostvtex) {
541:           Vr   = xarr[offset];
542:           Vi   = xarr[offset+1];

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

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

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

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

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

569:           e = connedges[i];
570:           DMNetworkGetComponentKeyOffset(networkdm,e,0,&keye,&offsetd);
571:           branch = (Branch*)(arr+offsetd);

573:           Yfti = branch->yft[0];
574:           Yftr = branch->yft[1];

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

578:           vfrom = cone[0];
579:           vto   = cone[1];

581:           DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
582:           DMNetworkGetVariableOffset(networkdm,vto,&offsetto);

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

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

608:           gen = (Gen*)(arr+offsetd);
609:           idx = offset + 2;

611:           /* Generator state variables */
612:           Eqp   = xarr[idx];
613:           Edp   = xarr[idx+1];
614:           delta = xarr[idx+2];
615:           w     = xarr[idx+3];
616:           Id    = xarr[idx+4];
617:           Iq    = xarr[idx+5];
618:           Efd   = xarr[idx+6];
619:           RF    = xarr[idx+7];
620:           VR    = xarr[idx+8];

622:           /* Generator parameters */
623:           Xd   = gen->Xd;
624:           Xdp  = gen->Xdp;
625:           Td0p = gen->Td0p;
626:           Xq   = gen->Xq;
627:           Xqp  = gen->Xqp;
628:           Tq0p = gen->Tq0p;
629:           TM   = gen->TM;
630:           D    = gen->D;
631:           M    = gen->M;
632:           Rs   = gen->Rs;
633:           k1   = gen->k1;
634:           k2   = gen->k2;
635:           KE   = gen->KE;
636:           TE   = gen->TE;
637:           TF   = gen->TF;
638:           KA   = gen->KA;
639:           KF   = gen->KF;
640:           Vref = gen->Vref;
641:           TA   = gen->TA;

643:           /* Generator differential equations */
644:           farr[idx]   = (Eqp + (Xd - Xdp)*Id - Efd)/Td0p + xdotarr[idx];
645:           farr[idx+1] = (Edp - (Xq - Xqp)*Iq)/Tq0p  + xdotarr[idx+1];
646:           farr[idx+2] = -w + W_S + xdotarr[idx+2];
647:           farr[idx+3] = (-TM + Edp*Id + Eqp*Iq + (Xqp - Xdp)*Id*Iq + D*(w - W_S))/M  + xdotarr[idx+3];

649:           Vr = xarr[offset]; /* Real part of generator terminal voltage */
650:           Vi = xarr[offset+1]; /* Imaginary part of the generator terminal voltage */

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

654:           /* Algebraic equations for stator currents */
655:           det = Rs*Rs + Xdp*Xqp;

657:           Zdq_inv[0] = Rs/det;
658:           Zdq_inv[1] = Xqp/det;
659:           Zdq_inv[2] = -Xdp/det;
660:           Zdq_inv[3] = Rs/det;

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

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

668:           farr[offset]   -= IGi;
669:           farr[offset+1] -= IGr;

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

674:           /* Exciter differential equations */
675:           farr[idx+6] = (KE*Efd + SE - VR)/TE + xdotarr[idx+6];
676:           farr[idx+7] = (RF - KF*Efd/TF)/TF + xdotarr[idx+7];
677:           farr[idx+8] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[idx+8];
678:         }
679:       } else if (key ==3){
680:         if (!ghostvtex) {
681:           PetscInt    k;
682:           PetscInt    ld_nsegsp;
683:           PetscInt    ld_nsegsq;
684:           PetscScalar *ld_alphap;
685:           PetscScalar *ld_betap,*ld_alphaq,*ld_betaq,PD0, QD0, IDr,IDi;

687:           load = (Load*)(arr+offsetd);

689:           /* Load Parameters */
690:           ld_nsegsp = load->ld_nsegsp;
691:           ld_alphap = load->ld_alphap;
692:           ld_betap  = load->ld_betap;
693:           ld_nsegsq = load->ld_nsegsq;
694:           ld_alphaq = load->ld_alphaq;
695:           ld_betaq  = load->ld_betaq;
696:           PD0       = load->PD0;
697:           QD0       = load->QD0;

699:           Vr  = xarr[offset]; /* Real part of generator terminal voltage */
700:           Vi  = xarr[offset+1]; /* Imaginary part of the generator terminal voltage */
701:           Vm  = PetscSqrtScalar(Vr*Vr + Vi*Vi);
702:           Vm2 = Vm*Vm;
703:           Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0);
704:           PD  = QD = 0.0;
705:           for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]);
706:           for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]);

708:           /* Load currents */
709:           IDr = (PD*Vr + QD*Vi)/Vm2;
710:           IDi = (-QD*Vr + PD*Vi)/Vm2;

712:           farr[offset]   += IDi;
713:           farr[offset+1] += IDr;
714:         }
715:       }
716:     }
717:   }

719:   VecRestoreArrayRead(localX,&xarr);
720:   VecRestoreArrayRead(localXdot,&xdotarr);
721:   VecRestoreArray(localF,&farr);
722:   DMRestoreLocalVector(networkdm,&localX);
723:   DMRestoreLocalVector(networkdm,&localXdot);

725:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
726:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
727:   DMRestoreLocalVector(networkdm,&localF);
728:   return(0);
729: }

731: /* This function is used for solving the algebraic system only during fault on and
732:    off times. It computes the entire F and then zeros out the part corresponding to
733:    differential equations
734:  F = [0;g(y)];
735: */
736: PetscErrorCode AlgFunction (SNES snes, Vec X, Vec F, void *ctx)
737: {
739:   DM             networkdm;
740:   Vec            localX,localF;
741:   PetscInt       vfrom,vto,offsetfrom,offsetto;
742:   PetscInt       v,vStart,vEnd,e;
743:   PetscScalar    *farr;
744:   Userctx        *user=(Userctx*)ctx;
745:   const PetscScalar *xarr;
746:   DMNetworkComponentGenericDataType *arr;

749:   VecSet(F,0.0);
750:   SNESGetDM(snes,&networkdm);
751:   DMGetLocalVector(networkdm,&localF);
752:   DMGetLocalVector(networkdm,&localX);
753:   VecSet(localF,0.0);

755:   /* update ghost values of locaX and locaXdot */
756:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
757:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

759:   VecGetArrayRead(localX,&xarr);
760:   VecGetArray(localF,&farr);

762:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
763:   DMNetworkGetComponentDataArray(networkdm,&arr);

765:   for (v=vStart; v < vEnd; v++) {
766:     PetscInt      i,j,offsetd,offset,key,numComps;
767:     PetscScalar   Vr, Vi, Yffr, Yffi, Vm, Vm2, Vm0, Vr0=0, Vi0=0, PD, QD;
768:     Bus           *bus;
769:     Gen           *gen;
770:     Load          *load;
771:     PetscBool     ghostvtex;

773:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
774:     DMNetworkGetNumComponents(networkdm,v,&numComps);
775:     DMNetworkGetVariableOffset(networkdm,v,&offset);

777:     for (j = 0; j < numComps; j++) {
778:       DMNetworkGetComponentKeyOffset(networkdm,v,j,&key,&offsetd);
779:       if (key == 1) {
780:         PetscInt       nconnedges;
781:         const PetscInt *connedges;

783:         bus = (Bus*)(arr+offsetd);
784:         if (!ghostvtex) {
785:           Vr = xarr[offset];
786:           Vi = xarr[offset+1];

788:           Yffr = bus->yff[1];
789:           Yffi = bus->yff[0];
790:           if (user->alg_flg){
791:             Yffr += user->ybusfault[bus->id*2+1];
792:             Yffi += user->ybusfault[bus->id*2];
793:           }
794:           Vr0 = bus->vr;
795:           Vi0 = bus->vi;

797:           farr[offset]   += Yffi*Vr + Yffr*Vi;
798:           farr[offset+1] += Yffr*Vr - Yffi*Vi;
799:         }
800:         DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);

802:         for (i=0; i < nconnedges; i++) {
803:           Branch         *branch;
804:           PetscInt       keye;
805:           PetscScalar    Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
806:           const PetscInt *cone;

808:           e = connedges[i];
809:           DMNetworkGetComponentKeyOffset(networkdm,e,0,&keye,&offsetd);
810:           branch = (Branch*)(arr+offsetd);

812:           Yfti = branch->yft[0];
813:           Yftr = branch->yft[1];

815:           DMNetworkGetConnectedVertices(networkdm,e,&cone);
816:           vfrom = cone[0];
817:           vto   = cone[1];

819:           DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
820:           DMNetworkGetVariableOffset(networkdm,vto,&offsetto);