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: /* 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:   PetscInt          i,j,row[1],col[2];
111:   PetscInt          *edgelist;
112:   PetscInt          nofgen[9] = {1,1,1,0,0,0,0,0,0}; /* Buses at which generators are incident */
113:   PetscInt          nofload[9] = {0,0,0,0,1,1,0,1,0}; /* Buses at which loads are incident */
114:   const PetscScalar  *varr;
115:   PetscScalar        M[3],D[3];
116:   Bus               *bus;
117:   Branch            *branch;
118:   Gen               *gen;
119:   Exc               *exc;
120:   Load              *load;
121:   Mat               Ybus;
122:   Vec               V0;

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

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

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

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

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

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

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

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

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

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

190:    VecGetArrayRead(V0,&varr);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

352:    VecRestoreArrayRead(V0,&varr);

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

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

380:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
381:   DMGetLocalVector(networkdm,&localX);

383:   VecSet(X,0.0);
384:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
385:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

387:   VecGetArray(localX,&xarr);

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

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

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

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

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

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

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

425:         Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta);
426:         Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta);

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

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

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

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

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

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

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

453:         xarr[offset] = Efd;
454:         xarr[offset+1] = RF;
455:         xarr[offset+2] = VR;

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

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

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

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

497:   VecSet(F,0.0);

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

654:           /* Add generator current injection to network */
655:           farr[offsetbus]   -= IGi;
656:           farr[offsetbus+1] -= IGr;

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

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

667:           Efd   = xarr[offsetexc];
668:           RF    = xarr[offsetexc+1];
669:           VR    = xarr[offsetexc+2];

671:           k1   = exc->k1;
672:           k2   = exc->k2;
673:           KE   = exc->KE;
674:           TE   = exc->TE;
675:           TF   = exc->TF;
676:           KA   = exc->KA;
677:           KF   = exc->KF;
678:           Vref = exc->Vref;
679:           TA   = exc->TA;

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

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

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

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

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

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

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

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

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

737:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
738:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
739:   DMRestoreLocalVector(networkdm,&localF);
740:   return 0;
741: }

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

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

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

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

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

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

783:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
784:     DMNetworkGetNumComponents(networkdm,v,&numComps);

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

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

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

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

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

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

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

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

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

line832">832: Vfr = xarr[offse
line833">833: Vfi = xarr[offsetf
line834">834: Vtr = xarr[off
line835">835: Vti = xarr[offse

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

line852">852: gen = (Gen*)(comp
line853">853: DMNetworkGetLocalVecOffset(networkdm,v,j,&offs

/* Generator state variables */
line856">856: Eqp = xarr[offs
line857">857: Edp = xarr[offset
line858">858: delta = xarr[offset
/* w = xarr[idx+3]; not being used */
line860">860: Id = xarr[offset
line861">861: Iq = xarr[offset

/* Generator parameters */
line864">864: Xdp = gen-&
line865">865: Xqp = gen-&
line866">866: Rs = gen-

/* Set generator differential equation residual functions to zero */
line869">869: farr[offsetgen]
line870">870: farr[offsetgen+
line871">871: farr[offsetgen+
line872">872: farr[offsetgen+

line874">874: ri2dq(Vr,Vi,delta,&Vd,&a

/* Algebraic equations for stator currents */
line877">877: det = Rs*Rs + X

line879">879: Zdq_inv[0] =
line880">880: Zdq_inv[1] = X
line881">881: Zdq_inv[2] = -X
line882">882: Zdq_inv[3] =

line884">884: farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq
line885">885: farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq

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

line890">890: farr[offsetbus]
line891">891: farr[offsetbus+1]

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

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

line910">910: load = (Load*)(comp

/* Load Parameters */
line913">913: ld_nsegsp = load->ld_
line914">914: ld_alphap = load->ld_
line915">915: ld_betap = load->ld
line916">916: ld_nsegsq = load->ld_
line917">917: ld_alphaq = load->ld_
line918">918: ld_betaq = load->ld

line920">920: PD0 = load-&
line921">921: QD0 = load-&

line923">923: Vm = PetscSqrtScalar(Vr*Vr +
line924">924: Vm2 =
line925">925: Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi
line926">926: PD = QD
line927">927: for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_bet
line928">928: for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_bet

/* Load currents */
line931">931: IDr = (PD*Vr + QD*V
line932">932: IDi = (-QD*Vr + PD*V

line934">934: farr[offsetbus]
line935">935: farr[offsetbus+1]
line936">936:
line937">937:
line938">938: line939">939:

line941">941: VecRestoreArrayRead(localX,&
line942">942: VecRestoreArray(localF,&
line943">943: DMRestoreLocalVector(networkdm,&l

line945">945: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES<
line946">946:
DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES<
line947">947:
DMRestoreLocalVector(networkdm,&l
line948">948: return line949">949

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

line975">975: PetscInitialize(&argc,&argv,"ex9busnetworkops"
line976">976: PetscOptionsGetInt(NULL,NULL,"-nc",&nc
line977">977: MPI_Comm_size(PETSC_COMM_WORLD,&
line978">978: MPI_Comm_rank(PETSC_COMM_WORLD,&

/* Read initial voltage vector and Ybus */
line981">981: if (rank
line982">982: read_data(nc,&gen,&exc,&load,&bus,&branch,&edg
line983">983:

line985">985: DMNetworkCreate(PETSC_COMM_WORLD,&netw
line986">986: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(Branch),&componentk
line987">987: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(Bus),&componentk
line988">988: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(Gen),&componentk
line989">989: DMNetworkRegisterComponent(networkdm,"excstruct",sizeof(Exc),&componentk
line990">990: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(Load),&componentk

line992">992: PetscLogStageRegister("Create network",&s
line993">993: PetscLogStagePush(s

/* Set local number of edges and edge connectivity */
line996">996: if (rank == 0) numEdges = NBRANCH*nc+
line997">997: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE<
line998">998:
DMNetworkAddSubnetwork(networkdm,NULL,numEdges,edgelist

/* Set up the network layout */
line1001">1001: DMNetworkLayoutSetUp(netw

line1003">1003: if (rank
line1004">1004: PetscFree(edg
line1005">1005:

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

line1017">1017: for (i = vStart; i < vEnd;
line1018">1018: DMNetworkAddComponent(networkdm,i,componentkey[1],&bus[i-vSta
line1019">1019: if (bus[i-vStart].no
line1020">1020: for (j = 0; j < bus[i-vStart].nofgen;
/* Add generator */
line1022">1022: DMNetworkAddComponent(networkdm,i,componentkey[2],&gen[genj
/* Add exciter */
line1024">1024: DMNetworkAddComponent(networkdm,i,componentkey[3],&exc[excj
line1025">1025:
line1026">1026:
line1027">1027: if (bus[i-vStart].nof
line1028">1028: for (j=0; j < bus[i-vStart].nofload;
line1029">1029: DMNetworkAddComponent(networkdm,i,componentkey[4],&load[loadj
line1030">1030:
line1031">1031:
line1032">1032: line1033">1033:

line1035">1035: DMSetUp(netw

line1037">1037: if (rank
line1038">1038: PetscFree5(bus,gen,load,branc
line1039">1039:

/* for parallel options: Network partitioning and distribution of data */
line1042">1042: if (size &g
line1043">1043: DMNetworkDistribute(&networ
line1044">1044:
line1045">1045: PetscLogStagePop

line1047">1047:
DMCreateGlobalVector(networkdm,&

line1049">1049: SetInitialGuess(networ

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

line1065">1065: PetscMalloc1(18*nc,&user.ybus
line1066">1066: for (i = 0; i < 18*nc;
line1067">1067: user.ybusfault[
line1068">1068:
line1069">1069: user.ybusfault[user.faultbus*2+1] = 1/user.
line1070">1070: PetscOptionsEnd

/* Setup
TS solver */
/*--------------------------------------------------------*/
line1074">1074: TSCreate(PETSC_COMM_WORLD,&a
line1075">1075: TSSetDM(ts,(DM)netw
line1076">1076: TSSetType(ts,TSC

line1078">1078:
TSGetSNES(ts,&
line1079">1079: SNESGetKSP(snes,&am
line1080">1080: KSPGetPC(ksp,&a
line1081">1081: PCSetType(pc,PCBJACOB

line1083">1083:
TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&
line1084">1084: TSSetMaxTime(ts,user.tfa
line1085">1085: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1086">1086:
TSSetTimeStep(ts
line1087">1087: 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. */
line1091">1091: user.alg_flg = PETSC_FAL

/* Prefault period */
line1094">1094:
if (rank
line1095">1095: PetscPrintf(PETSC_COMM_SELF,"... (1) Prefault period ... \n" line1096">1096:

line1098">1098: TSSetSolution
line1099">1099: TSSetUp line1100">1100: TSSolve

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

line1105">1105: SNESCreate(PETSC_COMM_WORLD,&sne
line1106">1106: SNESSetDM(snes_alg,(DM)netw
line1107">1107: SNESSetFunction(snes_alg,F_alg,AlgFunction,&
line1108">1108: SNESSetOptionsPrefix(snes_alg,"alg_" line1109">1109: SNESSetFromOptions(sne

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

/* Solve the algebraic equations */
line1117">1117:
if (rank
line1118">1118: PetscPrintf(PETSC_COMM_SELF,"\n... (2) Apply disturbance, solve algebraic equations ... \n" line1119">1119:
line1120">1120: SNESSolve(snes_alg,N

/* Disturbance period */
line1123">1123: TSSetTime(ts,user.tfa
line1124">1124: TSSetMaxTime(ts,user.tfau
line1125">1125: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1126">1126:
TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&

line1128">1128: user.alg_flg = PETSC_TR
line1129">1129:
if (rank
line1130">1130: PetscPrintf(PETSC_COMM_SELF,"\n... (3) Disturbance period ... \n" line1131">1131:
line1132">1132: TSSolve

/* Remove the fault */
line1135">1135: SNESSetFunction(snes_alg,F_alg,AlgFunction,&

line1137">1137: user.alg_flg = PETSC_FAL
/* Solve the algebraic equations */
line1139">1139:
if (rank
line1140">1140: PetscPrintf(PETSC_COMM_SELF,"\n... (4) Remove fault, solve algebraic equations ... \n" line1141">1141:
line1142">1142: SNESSolve(snes_alg,N
line1143">1143: SNESDestroy(&sne

/* Post-disturbance period */
line1146">1146: TSSetTime(ts,user.tfau
line1147">1147: TSSetMaxTime(ts,user
line1148">1148: TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVE
line1149">1149:
TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&

line1151">1151: user.alg_flg = PETSC_FAL
line1152">1152:
if (rank
line1153">1153: PetscPrintf(PETSC_COMM_SELF,"\n... (5) Post-disturbance period ... \n" line1154">1154:
line1155">1155: TSSolve

line1157">1157: PetscFree(user.ybus
line1158">1158: VecDestroy(&
line1159">1159: VecDestroy(&
line1160">1160: DMDestroy(&netw
line1161">1161: TSDestroy(&a
line1162">1162: PetscFinalize
line1163">1163:
return line1164">1164:

/*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*/