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:    /* Alocate 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;

655:         }
656:       } else if (key == 3) {
657:         if (!ghostvtex) {
658:           PetscScalar    k1,k2,KE,TE,TF,KA,KF,Vref,TA; /* Generator parameters */
659:           PetscScalar    Efd,RF,VR; /* Exciter variables */

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

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

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

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

681:           /* Exciter differential equations */
682:           farr[offsetexc] = (KE*Efd + SE - VR)/TE + xdotarr[offsetexc];
683:           farr[offsetexc+1] = (RF - KF*Efd/TF)/TF + xdotarr[offsetexc+1];
684:           farr[offsetexc+2] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[offsetexc+2];

686:         }
687:       } else if (key ==4) {
688:         if (!ghostvtex) {
689:           PetscInt    k;
690:           PetscInt    ld_nsegsp;
691:           PetscInt    ld_nsegsq;
692:           PetscScalar *ld_alphap;
693:           PetscScalar *ld_betap,*ld_alphaq,*ld_betaq,PD0, QD0, IDr,IDi;

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

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

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

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

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

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

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

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

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

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

767:   VecGetArrayRead(localX,&xarr);
768:   VecGetArray(localF,&farr);

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

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

780:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
781:     DMNetworkGetNumComponents(networkdm,v,&numComps);

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

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

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

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

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

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

818:           Yfti = branch->yft[0];
819:           Yftr = branch->yft[1];

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

825:           DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom);
826:           DMNetworkGetLocalVecOffset(networkdm,vto,0,&offsetto);

line829">829: Vfr = xarr[offse
line830">830: Vfi = xarr[offsetf
line831">831: Vtr = xarr[off
line832">832: Vti = xarr[offse

line834">834: if (vfrom
line835">835: farr[offsetfrom] += Yftr*Vti + Yf
line836">836: farr[offsetfrom+1] += Yftr*Vtr - Yf
line837">837: } else line838">838: farr[offsetto] += Yftr*Vfi + Yf
line839">839: farr[offsetto+1] += Yftr*Vfr - Yf
line840">840:
line841">841:
line842">842: } else if (key
line843">843: if (!ghost
/* Generator variables */
/* Generator dq axis currents */
line846">846: PetscScalar Vd,Vq,IGr,IGi,Zdq_inv[
/* Generator parameters */

line849">849: gen = (Gen*)(comp
line850">850: DMNetworkGetLocalVecOffset(networkdm,v,j,&offs

/* Generator state variables */
line853">853: Eqp = xarr[offs
line854">854: Edp = xarr[offset
line855">855: delta = xarr[offset
/* w = xarr[idx+3]; not being used */
line857">857: Id = xarr[offset
line858">858: Iq = xarr[offset

/* Generator parameters */
line861">861: Xdp = gen-&
line862">862: Xqp = gen-&
line863">863: Rs = gen-

/* Set generator differential equation residual functions to zero */
line866">866: farr[offsetgen]
line867">867: farr[offsetgen+
line868">868: farr[offsetgen+
line869">869: farr[offsetgen+

line871">871: ri2dq(Vr,Vi,delta,&Vd,&a

/* Algebraic equations for stator currents */
line874">874: det = Rs*Rs + X

line876">876: Zdq_inv[0] =
line877">877: Zdq_inv[1] = X
line878">878: Zdq_inv[2] = -X
line879">879: Zdq_inv[3] =

line881">881: farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq
line882">882: farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq

/* Add generator current injection to network */
line885">885: dq2ri(Id,Iq,delta,&IGr,&am

line887">887: farr[offsetbus]
line888">888: farr[offsetbus+1]

/* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/
/* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */

line892">892:
line893">893: } else if (key
line894">894: if (!ghost
line895">895: PetscInt off
line896">896: DMNetworkGetLocalVecOffset(networkdm,v,j,&offs
/* Set exciter differential equation residual functions equal to zero*/
line898">898: farr[offsetex
line899">899: farr[offsetexc+
line900">900: farr[offsetexc+
line901">901:
line902">902: } else if (key
line903">903: if (!ghost
line904">904: PetscInt k,ld_nsegsp,ld_
line905">905: PetscScalar *ld_alphap,*ld_betap,*ld_alphaq,*ld_betaq,PD0,QD0,I

line907">907: load = (Load*)(comp

/* Load Parameters */
line910">910: ld_nsegsp = load->ld_
line911">911: ld_alphap = load->ld_
line912">912: ld_betap = load->ld
line913">913: ld_nsegsq = load->ld_
line914">914: ld_alphaq = load->ld_
line915">915: ld_betaq = load->ld

line917">917: PD0 = load-&
line918">918: QD0 = load-&

line920">920: Vm = PetscSqrtScalar(Vr*Vr +
line921">921: Vm2 =
line922">922: Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi
line923">923: PD = QD
line924">924: for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_bet
line925">925: for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_bet

/* Load currents */
line928">928: IDr = (PD*Vr + QD*V
line929">929: IDi = (-QD*Vr + PD*V

line931">931: farr[offsetbus]
line932">932: farr[offsetbus+1]
line933">933:
line934">934:
line935">935: line936">936:

line938">938: VecRestoreArrayRead(localX,&
line939">939: VecRestoreArray(localF,&
line940">940: DMRestoreLocalVector(networkdm,&l

line942">942: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES<
line943">943:
DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES<
line944">944:
DMRestoreLocalVector(networkdm,&l
line945">945: return line946">946

line948">948: int main(int argc,char ** argv) line949">949
line951">951: PetscInt i,j,*edgelist= NULL,eStart,eEnd,vStar
line952">952: PetscInt genj,excj,loadj,component
/* No. of copies (default = 1) */
line954">954: PetscMPIInt siz
line955">955: Vec X
line956">956: TS
line957">957: SNES snes_al
line958">958: Bus
line959">959: Branch *
line960">960: Gen
line961">961: Exc
line962">962: Load
line963">963: DM net
line964">964: #if defined(PETSC_USE_LOG)
line965">965: PetscLogStage
line966">966: #endif
line967">967: Userctx
line968">968: KSP
line969">969: PC
line970">970: PetscInt numEdg

line972">972: PetscInitialize(&argc,&argv,"ex9busnetworkops"
line973">973: PetscOptionsGetInt(NULL,NULL,"-nc",&nc
line974">974: MPI_Comm_size(PETSC_COMM_WORLD,&
line975">975: MPI_Comm_rank(PETSC_COMM_WORLD,&

/* Read initial voltage vector and Ybus */
line978">978: if (rank
line979">979: read_data(nc,&gen,&exc,&load,&bus,&branch,&edg
line980">980:

line982">982: DMNetworkCreate(PETSC_COMM_WORLD,&netw
line983">983: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(Branch),&componentk
line984">984: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(Bus),&componentk
line985">985: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(Gen),&componentk
line986">986: DMNetworkRegisterComponent(networkdm,"excstruct",sizeof(Exc),&componentk
line987">987: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(Load),&componentk

line989">989: PetscLogStageRegister("Create network",&s
line990">990: PetscLogStagePush(s

/* Set local number of edges and edge connectivity */
line993">993: if (rank == 0) numEdges = NBRANCH*nc+
line994">994: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE<
line995">995:
DMNetworkAddSubnetwork(networkdm,NULL,numEdges,edgelist

/* Set up the network layout */
line998">998: DMNetworkLayoutSetUp(netw

line1000">1000: if (rank
line1001">1001: PetscFree(edg
line1002">1002:

/* Add network components (physical parameters of nodes and branches) and number of variables */
line1005">1005: if (rank
line1006">1006: DMNetworkGetEdgeRange(networkdm,&eStart,&
line1007">1007: genj=0; loadj=0;
line1008">1008: for (i = eStart; i < eEnd;
line1009">1009: DMNetworkAddComponent(networkdm,i,componentkey[0],&branch[i-eSta
line1010">1010:
line1012">1012: DMNetworkGetVertexRange(networkdm,&vStart,&

line1014">1014: for (i = vStart; i < vEnd;
line1015">1015: DMNetworkAddComponent(networkdm,i,componentkey[1],&bus[i-vSta
line1016">1016: if (bus[i-vStart].no
line1017">1017: for (j = 0; j < bus[i-vStart].nofgen;
/* Add generator */
line1019">1019: DMNetworkAddComponent(networkdm,i,componentkey[2],&gen[genj
/* Add exciter */
line1021">1021: DMNetworkAddComponent(networkdm,i,componentkey[3],&exc[excj
line1022">1022:
line1023">1023:
line1024">1024: if (bus[i-vStart].nof
line1025">1025: for (j=0; j < bus[i-vStart].nofload;
line1026">1026: DMNetworkAddComponent(networkdm,i,componentkey[4],&load[loadj
line1027">1027:
line1028">1028:
line1029">1029: line1030">1030:

line1032">1032: DMSetUp(netw

line1034">1034: if (rank
line1035">1035: PetscFree5(bus,gen,load,branc
line1036">1036:

/* for parallel options: Network partitioning and distribution of data */
line1039">1039: if (size &g
line1040">1040: DMNetworkDistribute(&networ
line1041">1041:
line1042">1042: PetscLogStagePop

line1044">1044:
DMCreateGlobalVector(networkdm,&

line1046">1046: SetInitialGuess(networ

/* Options for fault simulation */
line1049">1049: PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","" line1050">1050: user.tfaulton
line1051">1051: user.tfaultoff
line1052">1052: user.Rfault =
line1053">1053: user.faultbu
line1054">1054: PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton
line1055">1055: PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff
line1056">1056: PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus
line1057">1057: user.t0
line1058">1058: user.tmax
line1059">1059: PetscOptionsReal("-t0","","",user.t0,&user.t0
line1060">1060: PetscOptionsReal("-tmax","","",user.tmax,&user.tmax

line1062">1062: PetscMalloc1(18*nc,&user.ybus
line1063">1063: for (i = 0; i < 18*nc;
line1064">1064: user.ybusfault[
line1065">1065:
line1066">1066: user.ybusfault[user.faultbus*2+1] = 1/user.
line1067">1067: PetscOptionsEnd

/* Setup
TS solver */
/*--------------------------------------------------------*/
line1071">1071: TSCreate(PETSC_COMM_WORLD,&a
line1072">1072: TSSetDM(ts,(DM)netw
line1073">1073: TSSetType(ts,TSC

line1075">1075:
TSGetSNES(ts,&
line1076">1076: SNESGetKSP(snes,&am
line1077">1077: KSPGetPC(ksp,&a
line1078">1078: PCSetType(pc,PCBJACOB

line1080">1080:
TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1081">1081: TSSetMaxTime(ts,user.tfa
line1082">1082: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1083">1083:
TSSetTimeStep(ts
line1084">1084: TSSetFromOptions
/*user.alg_flg =
PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix.
eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */
line1088">1088: user.alg_flg = PETSC_FAL

/* Prefault period */
line1091">1091:
if (rank
line1092">1092: PetscPrintf(PETSC_COMM_SELF,"... (1) Prefault period ... \n" line1093">1093:

line1095">1095: TSSetSolution
line1096">1096: TSSetUp line1097">1097: TSSolve

/* Create the nonlinear solver for solving the algebraic system */
line1100">1100: VecDuplicate(X,&

line1102">1102: SNESCreate(PETSC_COMM_WORLD,&sne
line1103">1103: SNESSetDM(snes_alg,(DM)netw
line1104">1104: SNESSetFunction(snes_alg,F_alg,AlgFunction,&
line1105">1105: SNESSetOptionsPrefix(snes_alg,"alg_" line1106">1106: SNESSetFromOptions(sne

/* Apply disturbance - resistive fault at user.faultbus */
/* This is done by adding shunt conductance to the diagonal location
in the Ybus matrix */
line1111">1111: user.alg_flg = PETSC_TR

/* Solve the algebraic equations */
line1114">1114:
if (rank
line1115">1115: PetscPrintf(PETSC_COMM_SELF,"\n... (2) Apply disturbance, solve algebraic equations ... \n" line1116">1116:
line1117">1117: SNESSolve(snes_alg,N

/* Disturbance period */
line1120">1120: TSSetTime(ts,user.tfa
line1121">1121: TSSetMaxTime(ts,user.tfau
line1122">1122: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1123">1123:
TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&

line1125">1125: user.alg_flg = PETSC_TR
line1126">1126:
if (rank
line1127">1127: PetscPrintf(PETSC_COMM_SELF,"\n... (3) Disturbance period ... \n" line1128">1128:
line1129">1129: TSSolve

/* Remove the fault */
line1132">1132: SNESSetFunction(snes_alg,F_alg,AlgFunction,&

line1134">1134: user.alg_flg = PETSC_FAL
/* Solve the algebraic equations */
line1136">1136:
if (rank
line1137">1137: PetscPrintf(PETSC_COMM_SELF,"\n... (4) Remove fault, solve algebraic equations ... \n" line1138">1138:
line1139">1139: SNESSolve(snes_alg,N
line1140">1140: SNESDestroy(&sne

/* Post-disturbance period */
line1143">1143: TSSetTime(ts,user.tfau
line1144">1144: TSSetMaxTime(ts,user
line1145">1145: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1146">1146:
TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&

line1148">1148: user.alg_flg = PETSC_FAL
line1149">1149:
if (rank
line1150">1150: PetscPrintf(PETSC_COMM_SELF,"\n... (5) Post-disturbance period ... \n" line1151">1151:
line1152">1152: TSSolve

line1154">1154: PetscFree(user.ybus
line1155">1155: VecDestroy(&
line1156">1156: VecDestroy(&
line1157">1157: DMDestroy(&netw
line1158">1158: TSDestroy(&a
line1159">1159: PetscFinalize
line1160">1160:
return line1161">1161:

/*TEST

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

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

test:
suffix: 2
nsize: 2
args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason
localrunfiles: X.bin Ybus.bin ex9busnetworkops

TEST*/