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

 10: /* T
 11:    Concepts: DMNetwork
 12:    Concepts: PETSc TS solver

 14:    This example was modified from ex9busdmnetwork.c.
 15: */

 17:  #include <petscts.h>
 18:  #include <petscdmnetwork.h>

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

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

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

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

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

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

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

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

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

107: /* Used to read data into the DMNetwork components */
108: PetscErrorCode read_data(PetscInt nc, Gen **pgen,Exc **pexc, Load **pload,Bus **pbus, Branch **pbranch, PetscInt **pedgelist)
109: {
110:   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:   const PetscScalar  *varr;
116:   PetscScalar        M[3],D[3];
117:   Bus               *bus;
118:   Branch            *branch;
119:   Gen               *gen;
120:   Exc               *exc;
121:   Load              *load;
122:   Mat               Ybus;
123:   Vec               V0;

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

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

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

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

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

169:    MatCreate(PETSC_COMM_SELF,&Ybus);
170:    MatSetType(Ybus,MATBAIJ);
171:    MatLoad(Ybus,Ybusview);

173:    /* Destroy unnecessary stuff */
174:    PetscViewerDestroy(&Xview);
175:    PetscViewerDestroy(&Ybusview);

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

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

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

191:    VecGetArrayRead(V0,&varr);

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

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

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

242:    /* read load data */
243:    for (i = 0; i<nc; i++){
244:      for (j = 0; j < NLOAD; j++) {
245:        load[i*3+j].id        = i*3+j;
246:        load[i*3+j].PD0       = PD0[j];
247:        load[i*3+j].QD0       = QD0[j];
248:        load[i*3+j].ld_nsegsp = ld_nsegsp[j];

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

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

258:        load[i*3+j].ld_betap[0] = ld_betap[0];
259:        load[i*3+j].ld_betap[1] = ld_betap[1];
260:        load[i*3+j].ld_betap[2] = ld_betap[2];
261:        load[i*3+j].ld_nsegsq   = ld_nsegsq[j];

263:        load[i*3+j].ld_betaq[0] = ld_betaq[0];
264:        load[i*3+j].ld_betaq[1] = ld_betaq[1];
265:        load[i*3+j].ld_betaq[2] = ld_betaq[2];
266:      }
267:    }
268:    PetscCalloc1(2*NBRANCH*nc+2*(nc-1),&edgelist);

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

316:    /* 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 */
317:     for (i = 1; i<nc; i++){
318:         edgelist[18*nc+2*(i-1)]   = 8+(i-1)*9;
319:         edgelist[18*nc+2*(i-1)+1] = 9*i;

321:         /* adding admittances to the off-diagonal elements */
322:         branch[9*nc+(i-1)].id     = 9*nc+(i-1);
323:         branch[9*nc+(i-1)].yft[0] = 17.3611;
324:         branch[9*nc+(i-1)].yft[1] = -0.0301407;

326:         /* subtracting admittances from the diagonal elements */
327:         bus[9*i-1].yff[0] -= 17.3611;
328:         bus[9*i-1].yff[1] -= -0.0301407;

330:         bus[9*i].yff[0]   -= 17.3611;
331:         bus[9*i].yff[1]   -= -0.0301407;
332:     }

334:     /* read branch data */
335:     for (i = 0; i<nc; i++){
336:       for (j = 0; j < NBRANCH; j++) {
337:         branch[i*9+j].id  = i*9+j;

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

346:    *pgen      = gen;
347:    *pexc      = exc;
348:    *pload     = load;
349:    *pbus      = bus;
350:    *pbranch   = branch;
351:    *pedgelist = edgelist;

353:    VecRestoreArrayRead(V0,&varr);

355:    /* Destroy unnecessary stuff */
356:    MatDestroy(&Ybus);
357:    VecDestroy(&V0);
358:    return(0);
359: }

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

383:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
384:   DMGetLocalVector(networkdm,&localX);

386:   VecSet(X,0.0);
387:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
388:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

390:   VecGetArray(localX,&xarr);

392:   for (v = vStart; v < vEnd; v++) {
393:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
394:     if (ghostvtex) continue;

396:     DMNetworkGetNumComponents(networkdm,v,&numComps);
397:     for (j=0; j < numComps; j++) {
398:       DMNetworkGetComponent(networkdm,v,j,&key,&component);
399:       if (key == 1) {
400:         bus = (Bus*)(component);

402:         DMNetworkGetComponentVariableOffset(networkdm,v,j,&offset);
403:         xarr[offset]   = bus->vr;
404:         xarr[offset+1] = bus->vi;

406:         Vr = bus->vr;
407:         Vi = bus->vi;
408:       } else if(key == 2) {
409:         gen = (Gen*)(component);
410:         DMNetworkGetComponentVariableOffset(networkdm,v,j,&offset);
411:         Vm  = PetscSqrtScalar(Vr*Vr + Vi*Vi);
412:         Vm2 = Vm*Vm;
413:         /* Real part of gen current */
414:         IGr = (Vr*gen->PG + Vi*gen->QG)/Vm2;                                     
415:         /* Imaginary part of gen current */
416:         IGi = (Vi*gen->PG - Vr*gen->QG)/Vm2;

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

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

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

428:         Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta);
429:         Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta);

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

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

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

439:         xarr[offset]   = Eqp;
440:         xarr[offset+1] = Edp;
441:         xarr[offset+2] = delta;
442:         xarr[offset+3] = W_S;
443:         xarr[offset+4] = Id;
444:         xarr[offset+5] = Iq;

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

448:       } else if(key == 3) {
449:         exc = (Exc*)(component);
450:         DMNetworkGetComponentVariableOffset(networkdm,v,j,&offset);
451:         
452:         SE  = exc->k1*PetscExpScalar(exc->k2*Efd);
453:         VR  = exc->KE*Efd + SE;
454:         RF  = exc->KF*Efd/exc->TF;

456:         xarr[offset] = Efd;
457:         xarr[offset+1] = RF;
458:         xarr[offset+2] = VR;

460:         exc->Vref = Vm + (VR/exc->KA);
461:       }
462:     }
463:   }
464:   VecRestoreArray(localX,&xarr);
465:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
466:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
467:   DMRestoreLocalVector(networkdm,&localX);
468:   return(0);
469:  }

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

480: /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */
481: PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd,PetscScalar *Fq)
482: {
484:   *Fd =  Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta);
485:   *Fq =  Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta);
486:   return(0);
487: }

489: /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */
490: PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx *user)
491: {
492:   PetscErrorCode                     ierr;
493:   DM                                 networkdm;
494:   Vec                                localX,localXdot,localF;
495:   PetscInt                           vfrom,vto,offsetfrom,offsetto;
496:   PetscInt                           v,vStart,vEnd,e;
497:   PetscScalar                        *farr;
498:   PetscScalar                        Vd,Vq,SE;
499:   const PetscScalar                  *xarr,*xdotarr;
500:   void*                              component;
501:   PetscScalar                        Vr=0, Vi=0;

504:   VecSet(F,0.0);

506:   TSGetDM(ts,&networkdm);
507:   DMGetLocalVector(networkdm,&localF);
508:   DMGetLocalVector(networkdm,&localX);
509:   DMGetLocalVector(networkdm,&localXdot);
510:   VecSet(localF,0.0);

512:   /* update ghost values of localX and localXdot */
513:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
514:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

516:   DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot);
517:   DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot);

519:   VecGetArrayRead(localX,&xarr);
520:   VecGetArrayRead(localXdot,&xdotarr);
521:   VecGetArray(localF,&farr);

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

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

538:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
539:     DMNetworkGetNumComponents(networkdm,v,&numComps);

541:     for (j = 0; j < numComps; j++) {
542:       DMNetworkGetComponent(networkdm,v,j,&key,&component);
543:       if (key == 1) {
544:         PetscInt       nconnedges;
545:         const PetscInt *connedges;

547:         bus = (Bus*)(component);
548:         DMNetworkGetComponentVariableOffset(networkdm,v,j,&offsetbus);
549:         if (!ghostvtex) {
550:           Vr   = xarr[offsetbus];
551:           Vi   = xarr[offsetbus+1];

553:           Yffr = bus->yff[1];
554:           Yffi = bus->yff[0];

556:           if (user->alg_flg){
557:             Yffr += user->ybusfault[bus->id*2+1];
558:             Yffi += user->ybusfault[bus->id*2];
559:           }
560:           Vr0 = bus->vr;
561:           Vi0 = bus->vi;

563:           /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here.
564:            The generator current injection, IG, and load current injection, ID are added later
565:            */
566:           farr[offsetbus] +=  Yffi*Vr + Yffr*Vi; /* imaginary current due to diagonal elements */
567:           farr[offsetbus+1] += Yffr*Vr - Yffi*Vi; /* real current due to diagonal elements */
568:         }

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

572:         for (i=0; i < nconnedges; i++) {
573:           Branch         *branch;
574:           PetscInt       keye;
575:           PetscScalar    Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
576:           const PetscInt *cone;

578:           e = connedges[i];
579:           DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch);

581:           Yfti = branch->yft[0];
582:           Yftr = branch->yft[1];

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

586:           vfrom = cone[0];
587:           vto   = cone[1];

589:           DMNetworkGetComponentVariableOffset(networkdm,vfrom,0,&offsetfrom);
590:           DMNetworkGetComponentVariableOffset(networkdm,vto,0,&offsetto);

592:           /* From bus and to bus real and imaginary voltages */
593:           Vfr     = xarr[offsetfrom];
594:           Vfi     = xarr[offsetfrom+1];
595:           Vtr          = xarr[offsetto];
596:           Vti     = xarr[offsetto+1];

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

614:           gen = (Gen*)(component);
615:           DMNetworkGetComponentVariableOffset(networkdm,v,j,&offsetgen);

617:           /* Generator state variables */
618:           Eqp   = xarr[offsetgen];
619:           Edp   = xarr[offsetgen+1];
620:           delta = xarr[offsetgen+2];
621:           w     = xarr[offsetgen+3];
622:           Id    = xarr[offsetgen+4];
623:           Iq    = xarr[offsetgen+5];

625:           /* Generator parameters */
626:           Xd   = gen->Xd;
627:           Xdp  = gen->Xdp;
628:           Td0p = gen->Td0p;
629:           Xq   = gen->Xq;
630:           Xqp  = gen->Xqp;
631:           Tq0p = gen->Tq0p;
632:           TM   = gen->TM;
633:           D    = gen->D;
634:           M    = gen->M;
635:           Rs   = gen->Rs;

637:           DMNetworkGetComponentVariableOffset(networkdm,v,2,&offsetexc);
638:           Efd = xarr[offsetexc];

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

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

648:           /* Algebraic equations for stator currents */
649:           det = Rs*Rs + Xdp*Xqp;

651:           Zdq_inv[0] = Rs/det;
652:           Zdq_inv[1] = Xqp/det;
653:           Zdq_inv[2] = -Xdp/det;
654:           Zdq_inv[3] = Rs/det;

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

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

661:           /* Add generator current injection to network */
662:           farr[offsetbus]   -= IGi;
663:           farr[offsetbus+1] -= IGr;

665:         }
666:       } else if (key == 3) {
667:         if(!ghostvtex) {
668:           PetscScalar    k1,k2,KE,TE,TF,KA,KF,Vref,TA; /* Generator parameters */
669:           PetscScalar    Efd,RF,VR; /* Exciter variables */
670:           
671:           exc = (Exc*)(component);
672:           DMNetworkGetComponentVariableOffset(networkdm,v,j,&offsetexc);

674:           Efd   = xarr[offsetexc];
675:           RF    = xarr[offsetexc+1];
676:           VR    = xarr[offsetexc+2];

678:           k1   = exc->k1;
679:           k2   = exc->k2;
680:           KE   = exc->KE;
681:           TE   = exc->TE;
682:           TF   = exc->TF;
683:           KA   = exc->KA;
684:           KF   = exc->KF;
685:           Vref = exc->Vref;
686:           TA   = exc->TA;

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

691:           /* Exciter differential equations */
692:           farr[offsetexc] = (KE*Efd + SE - VR)/TE + xdotarr[offsetexc];
693:           farr[offsetexc+1] = (RF - KF*Efd/TF)/TF + xdotarr[offsetexc+1];
694:           farr[offsetexc+2] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[offsetexc+2];
695:           
696:         }
697:       } else if (key ==4){
698:         if (!ghostvtex) {
699:           PetscInt    k;
700:           PetscInt    ld_nsegsp;
701:           PetscInt    ld_nsegsq;
702:           PetscScalar *ld_alphap;
703:           PetscScalar *ld_betap,*ld_alphaq,*ld_betaq,PD0, QD0, IDr,IDi;

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

707:           /* Load Parameters */
708:           ld_nsegsp = load->ld_nsegsp;
709:           ld_alphap = load->ld_alphap;
710:           ld_betap  = load->ld_betap;
711:           ld_nsegsq = load->ld_nsegsq;
712:           ld_alphaq = load->ld_alphaq;
713:           ld_betaq  = load->ld_betaq;
714:           PD0       = load->PD0;
715:           QD0       = load->QD0;

717:           Vr  = xarr[offsetbus]; /* Real part of generator terminal voltage */
718:           Vi  = xarr[offsetbus+1]; /* Imaginary part of the generator terminal voltage */
719:           Vm  = PetscSqrtScalar(Vr*Vr + Vi*Vi); 
720:           Vm2 = Vm*Vm;
721:           Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0);
722:           PD  = QD = 0.0;
723:           for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]);
724:           for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]);

726:           /* Load currents */
727:           IDr = (PD*Vr + QD*Vi)/Vm2;
728:           IDi = (-QD*Vr + PD*Vi)/Vm2;

730:           /* Load current contribution to the network */
731:           farr[offsetbus]   += IDi;
732:           farr[offsetbus+1] += IDr;
733:         }
734:       }
735:     }
736:   }

738:   VecRestoreArrayRead(localX,&xarr);
739:   VecRestoreArrayRead(localXdot,&xdotarr);
740:   VecRestoreArray(localF,&farr);
741:   DMRestoreLocalVector(networkdm,&localX);
742:   DMRestoreLocalVector(networkdm,&localXdot);

744:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
745:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
746:   DMRestoreLocalVector(networkdm,&localF);
747:   return(0);
748: }

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

769:   VecSet(F,0.0);
770:   SNESGetDM(snes,&networkdm);
771:   DMGetLocalVector(networkdm,&localF);
772:   DMGetLocalVector(networkdm,&localX);
773:   VecSet(localF,0.0);

775:   /* update ghost values of locaX and locaXdot */
776:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
777:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

779:   VecGetArrayRead(localX,&xarr);
780:   VecGetArray(localF,&farr);

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

784:   for (v=vStart; v < vEnd; v++) {
785:     PetscInt      i,j,offsetbus,offsetgen,key,numComps;
786:     PetscScalar   Yffr, Yffi, Vm, Vm2, Vm0, Vr0=0, Vi0=0, PD, QD;
787:     Bus           *bus;
788:     Gen           *gen;
789:     Load          *load;
790:     PetscBool     ghostvtex;

792:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
793:     DMNetworkGetNumComponents(networkdm,v,&numComps);

795:     for (j = 0; j < numComps; j++) {
796:       DMNetworkGetComponent(networkdm,v,j,&key,&component);
797:       if (key == 1) {
798:         PetscInt       nconnedges;
799:         const PetscInt *connedges;

801:         bus = (Bus*)(component);
802:         DMNetworkGetComponentVariableOffset(networkdm,v,j,&offsetbus);
803:         if (!ghostvtex) {
804:           Vr = xarr[offsetbus];
805:           Vi = xarr[offsetbus+1];

807:           Yffr = bus->yff[1];
808:           Yffi = bus->yff[0];
809:           if (user->alg_flg){
810:             Yffr += user->ybusfault[bus->id*2+1];
811:             Yffi += user->ybusfault[bus->id*2];
812:           }
813:           Vr0 = bus->vr;
814:           Vi0 = bus->vi;

816:           farr[offsetbus]   += Yffi*Vr + Yffr*Vi;
817:           farr[offsetbus+1] += Yffr*Vr - Yffi*Vi;
818:         }
819:         DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);

821:         for (i=0; i < nconnedges; i++) {
822:           Branch         *branch;
823:           PetscInt       keye;
824:           PetscScalar    Yfti, Yftr, Vfr, Vfi, Vtr, Vti;
825:           const PetscInt *cone;

827:           e = connedges[i];
828:           DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch);

830:           Yfti = branch->yft[0];
831:           Yftr = branch->yft[1];

833:           DMNetworkGetConnectedVertices(networkdm,e,&cone);
834:           vfrom = cone[0];
835:           vto   = cone[1];

837:           DMNetworkGetComponentVariableOffset(networkdm,vfrom,0,&offsetfrom);
838:           DMNetworkGetComponentVariableOffset(networkdm,vto,0,&offsetto);