1 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"; 9 10 /* T 11 Concepts: DMNetwork 12 Concepts: PETSc TS solver 13 14 This example was modified from ex9busdmnetwork.c. 15 */ 16 17 #include <petscts.h> 18 #include <petscdmnetwork.h> 19 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 */ 26 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 */ 32 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; 46 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; 58 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; 67 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 81 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; 91 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; 96 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; 106 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; 123 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}; 128 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 */ 138 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) */ 148 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; 160 161 PetscFunctionBeginUser; 162 /* Read V0 and Ybus from files */ 163 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF,"X.bin",FILE_MODE_READ,&Xview)); 164 PetscCall(PetscViewerBinaryOpen(PETSC_COMM_SELF,"Ybus.bin",FILE_MODE_READ,&Ybusview)); 165 PetscCall(VecCreate(PETSC_COMM_SELF,&V0)); 166 PetscCall(VecLoad(V0,Xview)); 167 168 PetscCall(MatCreate(PETSC_COMM_SELF,&Ybus)); 169 PetscCall(MatSetType(Ybus,MATBAIJ)); 170 PetscCall(MatLoad(Ybus,Ybusview)); 171 172 /* Destroy unnecessary stuff */ 173 PetscCall(PetscViewerDestroy(&Xview)); 174 PetscCall(PetscViewerDestroy(&Ybusview)); 175 176 PetscCall(MatGetLocalSize(Ybus,&m,&n)); 177 neqs_net = 2*NBUS; /* # eqs. for network subsystem */ 178 PetscCheckFalse(m != neqs_net || n != neqs_net,PETSC_COMM_SELF,PETSC_ERR_ARG_SIZ,"matrix Ybus is in wrong sizes"); 179 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]; 186 187 /* Alocate memory for bus, generators, exciter, loads and branches */ 188 PetscCall(PetscCalloc5(NBUS*nc,&bus,NGEN*nc,&gen,NLOAD*nc,&load,NBRANCH*nc+(nc-1),&branch,NGEN*nc,&exc)); 189 190 PetscCall(VecGetArrayRead(V0,&varr)); 191 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 PetscCall(MatGetValues(Ybus,1,row,2,col,bus[i*9+j].yff)); 205 } 206 } 207 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 } 226 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 } 240 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]; 248 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]; 252 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]; 256 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]; 261 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 PetscCall(PetscCalloc1(2*NBRANCH*nc+2*(nc-1),&edgelist)); 268 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 } 314 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; 319 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; 324 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; 328 329 bus[9*i].yff[0] -= 17.3611; 330 bus[9*i].yff[1] -= -0.0301407; 331 } 332 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; 337 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 PetscCall(MatGetValues(Ybus,1,row,2,col,branch[i*9+j].yft));/*imaginary part of admittance*/ 342 } 343 } 344 345 *pgen = gen; 346 *pexc = exc; 347 *pload = load; 348 *pbus = bus; 349 *pbranch = branch; 350 *pedgelist = edgelist; 351 352 PetscCall(VecRestoreArrayRead(V0,&varr)); 353 354 /* Destroy unnecessary stuff */ 355 PetscCall(MatDestroy(&Ybus)); 356 PetscCall(VecDestroy(&V0)); 357 PetscFunctionReturn(0); 358 } 359 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; 379 380 PetscFunctionBegin; 381 PetscCall(DMNetworkGetVertexRange(networkdm,&vStart,&vEnd)); 382 PetscCall(DMGetLocalVector(networkdm,&localX)); 383 384 PetscCall(VecSet(X,0.0)); 385 PetscCall(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 386 PetscCall(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 387 388 PetscCall(VecGetArray(localX,&xarr)); 389 390 for (v = vStart; v < vEnd; v++) { 391 PetscCall(DMNetworkIsGhostVertex(networkdm,v,&ghostvtex)); 392 if (ghostvtex) continue; 393 394 PetscCall(DMNetworkGetNumComponents(networkdm,v,&numComps)); 395 for (j=0; j < numComps; j++) { 396 PetscCall(DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL)); 397 if (key == 1) { 398 bus = (Bus*)(component); 399 400 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offset)); 401 xarr[offset] = bus->vr; 402 xarr[offset+1] = bus->vi; 403 404 Vr = bus->vr; 405 Vi = bus->vi; 406 } else if (key == 2) { 407 gen = (Gen*)(component); 408 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offset)); 409 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 410 Vm2 = Vm*Vm; 411 /* Real part of gen current */ 412 IGr = (Vr*gen->PG + Vi*gen->QG)/Vm2; 413 /* Imaginary part of gen current */ 414 IGi = (Vi*gen->PG - Vr*gen->QG)/Vm2; 415 416 /* Machine angle */ 417 delta = atan2(Vi+gen->Xq*IGr,Vr-gen->Xq*IGi); 418 theta = PETSC_PI/2.0 - delta; 419 420 /* d-axis stator current */ 421 Id = IGr*PetscCosScalar(theta) - IGi*PetscSinScalar(theta); 422 423 /* q-axis stator current */ 424 Iq = IGr*PetscSinScalar(theta) + IGi*PetscCosScalar(theta); 425 426 Vd = Vr*PetscCosScalar(theta) - Vi*PetscSinScalar(theta); 427 Vq = Vr*PetscSinScalar(theta) + Vi*PetscCosScalar(theta); 428 429 /* d-axis transient EMF */ 430 Edp = Vd + gen->Rs*Id - gen->Xqp*Iq; 431 432 /* q-axis transient EMF */ 433 Eqp = Vq + gen->Rs*Iq + gen->Xdp*Id; 434 435 gen->TM = gen->PG; 436 437 xarr[offset] = Eqp; 438 xarr[offset+1] = Edp; 439 xarr[offset+2] = delta; 440 xarr[offset+3] = W_S; 441 xarr[offset+4] = Id; 442 xarr[offset+5] = Iq; 443 444 Efd = Eqp + (gen->Xd - gen->Xdp)*Id; 445 446 } else if (key == 3) { 447 exc = (Exc*)(component); 448 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offset)); 449 450 SE = exc->k1*PetscExpScalar(exc->k2*Efd); 451 VR = exc->KE*Efd + SE; 452 RF = exc->KF*Efd/exc->TF; 453 454 xarr[offset] = Efd; 455 xarr[offset+1] = RF; 456 xarr[offset+2] = VR; 457 458 exc->Vref = Vm + (VR/exc->KA); 459 } 460 } 461 } 462 PetscCall(VecRestoreArray(localX,&xarr)); 463 PetscCall(DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X)); 464 PetscCall(DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X)); 465 PetscCall(DMRestoreLocalVector(networkdm,&localX)); 466 PetscFunctionReturn(0); 467 } 468 469 /* Converts from machine frame (dq) to network (phase a real,imag) reference frame */ 470 PetscErrorCode dq2ri(PetscScalar Fd,PetscScalar Fq,PetscScalar delta,PetscScalar *Fr,PetscScalar *Fi) 471 { 472 PetscFunctionBegin; 473 *Fr = Fd*PetscSinScalar(delta) + Fq*PetscCosScalar(delta); 474 *Fi = -Fd*PetscCosScalar(delta) + Fq*PetscSinScalar(delta); 475 PetscFunctionReturn(0); 476 } 477 478 /* Converts from network frame ([phase a real,imag) to machine (dq) reference frame */ 479 PetscErrorCode ri2dq(PetscScalar Fr,PetscScalar Fi,PetscScalar delta,PetscScalar *Fd,PetscScalar *Fq) 480 { 481 PetscFunctionBegin; 482 *Fd = Fr*PetscSinScalar(delta) - Fi*PetscCosScalar(delta); 483 *Fq = Fr*PetscCosScalar(delta) + Fi*PetscSinScalar(delta); 484 PetscFunctionReturn(0); 485 } 486 487 /* Computes F(t,U,U_t) where F() = 0 is the DAE to be solved. */ 488 PetscErrorCode FormIFunction(TS ts,PetscReal t,Vec X,Vec Xdot,Vec F,Userctx *user) 489 { 490 DM networkdm; 491 Vec localX,localXdot,localF; 492 PetscInt vfrom,vto,offsetfrom,offsetto; 493 PetscInt v,vStart,vEnd,e; 494 PetscScalar *farr; 495 PetscScalar Vd=0,Vq=0,SE; 496 const PetscScalar *xarr,*xdotarr; 497 void* component; 498 PetscScalar Vr=0, Vi=0; 499 500 PetscFunctionBegin; 501 PetscCall(VecSet(F,0.0)); 502 503 PetscCall(TSGetDM(ts,&networkdm)); 504 PetscCall(DMGetLocalVector(networkdm,&localF)); 505 PetscCall(DMGetLocalVector(networkdm,&localX)); 506 PetscCall(DMGetLocalVector(networkdm,&localXdot)); 507 PetscCall(VecSet(localF,0.0)); 508 509 /* update ghost values of localX and localXdot */ 510 PetscCall(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 511 PetscCall(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 512 513 PetscCall(DMGlobalToLocalBegin(networkdm,Xdot,INSERT_VALUES,localXdot)); 514 PetscCall(DMGlobalToLocalEnd(networkdm,Xdot,INSERT_VALUES,localXdot)); 515 516 PetscCall(VecGetArrayRead(localX,&xarr)); 517 PetscCall(VecGetArrayRead(localXdot,&xdotarr)); 518 PetscCall(VecGetArray(localF,&farr)); 519 520 PetscCall(DMNetworkGetVertexRange(networkdm,&vStart,&vEnd)); 521 522 for (v=vStart; v < vEnd; v++) { 523 PetscInt i,j,offsetbus,offsetgen,offsetexc,key; 524 Bus *bus; 525 Gen *gen; 526 Exc *exc; 527 Load *load; 528 PetscBool ghostvtex; 529 PetscInt numComps; 530 PetscScalar Yffr,Yffi; /* Real and imaginary fault admittances */ 531 PetscScalar Vm,Vm2,Vm0; 532 PetscScalar Vr0=0,Vi0=0; 533 PetscScalar PD,QD; 534 535 PetscCall(DMNetworkIsGhostVertex(networkdm,v,&ghostvtex)); 536 PetscCall(DMNetworkGetNumComponents(networkdm,v,&numComps)); 537 538 for (j = 0; j < numComps; j++) { 539 PetscCall(DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL)); 540 if (key == 1) { 541 PetscInt nconnedges; 542 const PetscInt *connedges; 543 544 bus = (Bus*)(component); 545 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetbus)); 546 if (!ghostvtex) { 547 Vr = xarr[offsetbus]; 548 Vi = xarr[offsetbus+1]; 549 550 Yffr = bus->yff[1]; 551 Yffi = bus->yff[0]; 552 553 if (user->alg_flg) { 554 Yffr += user->ybusfault[bus->id*2+1]; 555 Yffi += user->ybusfault[bus->id*2]; 556 } 557 Vr0 = bus->vr; 558 Vi0 = bus->vi; 559 560 /* Network current balance residual IG + Y*V + IL = 0. Only YV is added here. 561 The generator current injection, IG, and load current injection, ID are added later 562 */ 563 farr[offsetbus] += Yffi*Vr + Yffr*Vi; /* imaginary current due to diagonal elements */ 564 farr[offsetbus+1] += Yffr*Vr - Yffi*Vi; /* real current due to diagonal elements */ 565 } 566 567 PetscCall(DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges)); 568 569 for (i=0; i < nconnedges; i++) { 570 Branch *branch; 571 PetscInt keye; 572 PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti; 573 const PetscInt *cone; 574 575 e = connedges[i]; 576 PetscCall(DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL)); 577 578 Yfti = branch->yft[0]; 579 Yftr = branch->yft[1]; 580 581 PetscCall(DMNetworkGetConnectedVertices(networkdm,e,&cone)); 582 583 vfrom = cone[0]; 584 vto = cone[1]; 585 586 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom)); 587 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vto,0,&offsetto)); 588 589 /* From bus and to bus real and imaginary voltages */ 590 Vfr = xarr[offsetfrom]; 591 Vfi = xarr[offsetfrom+1]; 592 Vtr = xarr[offsetto]; 593 Vti = xarr[offsetto+1]; 594 595 if (vfrom == v) { 596 farr[offsetfrom] += Yftr*Vti + Yfti*Vtr; 597 farr[offsetfrom+1] += Yftr*Vtr - Yfti*Vti; 598 } else { 599 farr[offsetto] += Yftr*Vfi + Yfti*Vfr; 600 farr[offsetto+1] += Yftr*Vfr - Yfti*Vfi; 601 } 602 } 603 } else if (key == 2) { 604 if (!ghostvtex) { 605 PetscScalar Eqp,Edp,delta,w; /* Generator variables */ 606 PetscScalar Efd; /* Exciter field voltage */ 607 PetscScalar Id,Iq; /* Generator dq axis currents */ 608 PetscScalar IGr,IGi,Zdq_inv[4],det; 609 PetscScalar Xd,Xdp,Td0p,Xq,Xqp,Tq0p,TM,D,M,Rs; /* Generator parameters */ 610 611 gen = (Gen*)(component); 612 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetgen)); 613 614 /* Generator state variables */ 615 Eqp = xarr[offsetgen]; 616 Edp = xarr[offsetgen+1]; 617 delta = xarr[offsetgen+2]; 618 w = xarr[offsetgen+3]; 619 Id = xarr[offsetgen+4]; 620 Iq = xarr[offsetgen+5]; 621 622 /* Generator parameters */ 623 Xd = gen->Xd; 624 Xdp = gen->Xdp; 625 Td0p = gen->Td0p; 626 Xq = gen->Xq; 627 Xqp = gen->Xqp; 628 Tq0p = gen->Tq0p; 629 TM = gen->TM; 630 D = gen->D; 631 M = gen->M; 632 Rs = gen->Rs; 633 634 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,2,&offsetexc)); 635 Efd = xarr[offsetexc]; 636 637 /* Generator differential equations */ 638 farr[offsetgen] = (Eqp + (Xd - Xdp)*Id - Efd)/Td0p + xdotarr[offsetgen]; 639 farr[offsetgen+1] = (Edp - (Xq - Xqp)*Iq)/Tq0p + xdotarr[offsetgen+1]; 640 farr[offsetgen+2] = -w + W_S + xdotarr[offsetgen+2]; 641 farr[offsetgen+3] = (-TM + Edp*Id + Eqp*Iq + (Xqp - Xdp)*Id*Iq + D*(w - W_S))/M + xdotarr[offsetgen+3]; 642 643 PetscCall(ri2dq(Vr,Vi,delta,&Vd,&Vq)); 644 645 /* Algebraic equations for stator currents */ 646 det = Rs*Rs + Xdp*Xqp; 647 648 Zdq_inv[0] = Rs/det; 649 Zdq_inv[1] = Xqp/det; 650 Zdq_inv[2] = -Xdp/det; 651 Zdq_inv[3] = Rs/det; 652 653 farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 654 farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 655 656 PetscCall(dq2ri(Id,Iq,delta,&IGr,&IGi)); 657 658 /* Add generator current injection to network */ 659 farr[offsetbus] -= IGi; 660 farr[offsetbus+1] -= IGr; 661 662 } 663 } else if (key == 3) { 664 if (!ghostvtex) { 665 PetscScalar k1,k2,KE,TE,TF,KA,KF,Vref,TA; /* Generator parameters */ 666 PetscScalar Efd,RF,VR; /* Exciter variables */ 667 668 exc = (Exc*)(component); 669 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetexc)); 670 671 Efd = xarr[offsetexc]; 672 RF = xarr[offsetexc+1]; 673 VR = xarr[offsetexc+2]; 674 675 k1 = exc->k1; 676 k2 = exc->k2; 677 KE = exc->KE; 678 TE = exc->TE; 679 TF = exc->TF; 680 KA = exc->KA; 681 KF = exc->KF; 682 Vref = exc->Vref; 683 TA = exc->TA; 684 685 Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq); 686 SE = k1*PetscExpScalar(k2*Efd); 687 688 /* Exciter differential equations */ 689 farr[offsetexc] = (KE*Efd + SE - VR)/TE + xdotarr[offsetexc]; 690 farr[offsetexc+1] = (RF - KF*Efd/TF)/TF + xdotarr[offsetexc+1]; 691 farr[offsetexc+2] = (VR - KA*RF + KA*KF*Efd/TF - KA*(Vref - Vm))/TA + xdotarr[offsetexc+2]; 692 693 } 694 } else if (key ==4) { 695 if (!ghostvtex) { 696 PetscInt k; 697 PetscInt ld_nsegsp; 698 PetscInt ld_nsegsq; 699 PetscScalar *ld_alphap; 700 PetscScalar *ld_betap,*ld_alphaq,*ld_betaq,PD0, QD0, IDr,IDi; 701 702 load = (Load*)(component); 703 704 /* Load Parameters */ 705 ld_nsegsp = load->ld_nsegsp; 706 ld_alphap = load->ld_alphap; 707 ld_betap = load->ld_betap; 708 ld_nsegsq = load->ld_nsegsq; 709 ld_alphaq = load->ld_alphaq; 710 ld_betaq = load->ld_betaq; 711 PD0 = load->PD0; 712 QD0 = load->QD0; 713 714 Vr = xarr[offsetbus]; /* Real part of generator terminal voltage */ 715 Vi = xarr[offsetbus+1]; /* Imaginary part of the generator terminal voltage */ 716 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 717 Vm2 = Vm*Vm; 718 Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0); 719 PD = QD = 0.0; 720 for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]); 721 for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 722 723 /* Load currents */ 724 IDr = (PD*Vr + QD*Vi)/Vm2; 725 IDi = (-QD*Vr + PD*Vi)/Vm2; 726 727 /* Load current contribution to the network */ 728 farr[offsetbus] += IDi; 729 farr[offsetbus+1] += IDr; 730 } 731 } 732 } 733 } 734 735 PetscCall(VecRestoreArrayRead(localX,&xarr)); 736 PetscCall(VecRestoreArrayRead(localXdot,&xdotarr)); 737 PetscCall(VecRestoreArray(localF,&farr)); 738 PetscCall(DMRestoreLocalVector(networkdm,&localX)); 739 PetscCall(DMRestoreLocalVector(networkdm,&localXdot)); 740 741 PetscCall(DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F)); 742 PetscCall(DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F)); 743 PetscCall(DMRestoreLocalVector(networkdm,&localF)); 744 PetscFunctionReturn(0); 745 } 746 747 /* This function is used for solving the algebraic system only during fault on and 748 off times. It computes the entire F and then zeros out the part corresponding to 749 differential equations 750 F = [0;g(y)]; 751 */ 752 PetscErrorCode AlgFunction (SNES snes, Vec X, Vec F, void *ctx) 753 { 754 DM networkdm; 755 Vec localX,localF; 756 PetscInt vfrom,vto,offsetfrom,offsetto; 757 PetscInt v,vStart,vEnd,e; 758 PetscScalar *farr; 759 Userctx *user=(Userctx*)ctx; 760 const PetscScalar *xarr; 761 void* component; 762 PetscScalar Vr=0,Vi=0; 763 764 PetscFunctionBegin; 765 PetscCall(VecSet(F,0.0)); 766 PetscCall(SNESGetDM(snes,&networkdm)); 767 PetscCall(DMGetLocalVector(networkdm,&localF)); 768 PetscCall(DMGetLocalVector(networkdm,&localX)); 769 PetscCall(VecSet(localF,0.0)); 770 771 /* update ghost values of locaX and locaXdot */ 772 PetscCall(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 773 PetscCall(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 774 775 PetscCall(VecGetArrayRead(localX,&xarr)); 776 PetscCall(VecGetArray(localF,&farr)); 777 778 PetscCall(DMNetworkGetVertexRange(networkdm,&vStart,&vEnd)); 779 780 for (v=vStart; v < vEnd; v++) { 781 PetscInt i,j,offsetbus,offsetgen,key,numComps; 782 PetscScalar Yffr, Yffi, Vm, Vm2, Vm0, Vr0=0, Vi0=0, PD, QD; 783 Bus *bus; 784 Gen *gen; 785 Load *load; 786 PetscBool ghostvtex; 787 788 PetscCall(DMNetworkIsGhostVertex(networkdm,v,&ghostvtex)); 789 PetscCall(DMNetworkGetNumComponents(networkdm,v,&numComps)); 790 791 for (j = 0; j < numComps; j++) { 792 PetscCall(DMNetworkGetComponent(networkdm,v,j,&key,&component,NULL)); 793 if (key == 1) { 794 PetscInt nconnedges; 795 const PetscInt *connedges; 796 797 bus = (Bus*)(component); 798 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetbus)); 799 if (!ghostvtex) { 800 Vr = xarr[offsetbus]; 801 Vi = xarr[offsetbus+1]; 802 803 Yffr = bus->yff[1]; 804 Yffi = bus->yff[0]; 805 if (user->alg_flg) { 806 Yffr += user->ybusfault[bus->id*2+1]; 807 Yffi += user->ybusfault[bus->id*2]; 808 } 809 Vr0 = bus->vr; 810 Vi0 = bus->vi; 811 812 farr[offsetbus] += Yffi*Vr + Yffr*Vi; 813 farr[offsetbus+1] += Yffr*Vr - Yffi*Vi; 814 } 815 PetscCall(DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges)); 816 817 for (i=0; i < nconnedges; i++) { 818 Branch *branch; 819 PetscInt keye; 820 PetscScalar Yfti, Yftr, Vfr, Vfi, Vtr, Vti; 821 const PetscInt *cone; 822 823 e = connedges[i]; 824 PetscCall(DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL)); 825 826 Yfti = branch->yft[0]; 827 Yftr = branch->yft[1]; 828 829 PetscCall(DMNetworkGetConnectedVertices(networkdm,e,&cone)); 830 vfrom = cone[0]; 831 vto = cone[1]; 832 833 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vfrom,0,&offsetfrom)); 834 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vto,0,&offsetto)); 835 836 /*From bus and to bus real and imaginary voltages */ 837 Vfr = xarr[offsetfrom]; 838 Vfi = xarr[offsetfrom+1]; 839 Vtr = xarr[offsetto]; 840 Vti = xarr[offsetto+1]; 841 842 if (vfrom == v) { 843 farr[offsetfrom] += Yftr*Vti + Yfti*Vtr; 844 farr[offsetfrom+1] += Yftr*Vtr - Yfti*Vti; 845 } else { 846 farr[offsetto] += Yftr*Vfi + Yfti*Vfr; 847 farr[offsetto+1] += Yftr*Vfr - Yfti*Vfi; 848 } 849 } 850 } else if (key == 2) { 851 if (!ghostvtex) { 852 PetscScalar Eqp,Edp,delta; /* Generator variables */ 853 PetscScalar Id,Iq; /* Generator dq axis currents */ 854 PetscScalar Vd,Vq,IGr,IGi,Zdq_inv[4],det; 855 PetscScalar Xdp,Xqp,Rs; /* Generator parameters */ 856 857 gen = (Gen*)(component); 858 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetgen)); 859 860 /* Generator state variables */ 861 Eqp = xarr[offsetgen]; 862 Edp = xarr[offsetgen+1]; 863 delta = xarr[offsetgen+2]; 864 /* w = xarr[idx+3]; not being used */ 865 Id = xarr[offsetgen+4]; 866 Iq = xarr[offsetgen+5]; 867 868 /* Generator parameters */ 869 Xdp = gen->Xdp; 870 Xqp = gen->Xqp; 871 Rs = gen->Rs; 872 873 /* Set generator differential equation residual functions to zero */ 874 farr[offsetgen] = 0; 875 farr[offsetgen+1] = 0; 876 farr[offsetgen+2] = 0; 877 farr[offsetgen+3] = 0; 878 879 PetscCall(ri2dq(Vr,Vi,delta,&Vd,&Vq)); 880 881 /* Algebraic equations for stator currents */ 882 det = Rs*Rs + Xdp*Xqp; 883 884 Zdq_inv[0] = Rs/det; 885 Zdq_inv[1] = Xqp/det; 886 Zdq_inv[2] = -Xdp/det; 887 Zdq_inv[3] = Rs/det; 888 889 farr[offsetgen+4] = Zdq_inv[0]*(-Edp + Vd) + Zdq_inv[1]*(-Eqp + Vq) + Id; 890 farr[offsetgen+5] = Zdq_inv[2]*(-Edp + Vd) + Zdq_inv[3]*(-Eqp + Vq) + Iq; 891 892 /* Add generator current injection to network */ 893 PetscCall(dq2ri(Id,Iq,delta,&IGr,&IGi)); 894 895 farr[offsetbus] -= IGi; 896 farr[offsetbus+1] -= IGr; 897 898 /* Vm = PetscSqrtScalar(Vd*Vd + Vq*Vq);*/ /* a compiler warning: "Value stored to 'Vm' is never read" - comment out by Hong Zhang */ 899 900 } 901 } else if (key == 3) { 902 if (!ghostvtex) { 903 PetscInt offsetexc; 904 PetscCall(DMNetworkGetLocalVecOffset(networkdm,v,j,&offsetexc)); 905 /* Set exciter differential equation residual functions equal to zero*/ 906 farr[offsetexc] = 0; 907 farr[offsetexc+1] = 0; 908 farr[offsetexc+2] = 0; 909 } 910 } else if (key == 4) { 911 if (!ghostvtex) { 912 PetscInt k,ld_nsegsp,ld_nsegsq; 913 PetscScalar *ld_alphap,*ld_betap,*ld_alphaq,*ld_betaq,PD0,QD0,IDr,IDi; 914 915 load = (Load*)(component); 916 917 /* Load Parameters */ 918 ld_nsegsp = load->ld_nsegsp; 919 ld_alphap = load->ld_alphap; 920 ld_betap = load->ld_betap; 921 ld_nsegsq = load->ld_nsegsq; 922 ld_alphaq = load->ld_alphaq; 923 ld_betaq = load->ld_betaq; 924 925 PD0 = load->PD0; 926 QD0 = load->QD0; 927 928 Vm = PetscSqrtScalar(Vr*Vr + Vi*Vi); 929 Vm2 = Vm*Vm; 930 Vm0 = PetscSqrtScalar(Vr0*Vr0 + Vi0*Vi0); 931 PD = QD = 0.0; 932 for (k=0; k < ld_nsegsp; k++) PD += ld_alphap[k]*PD0*PetscPowScalar((Vm/Vm0),ld_betap[k]); 933 for (k=0; k < ld_nsegsq; k++) QD += ld_alphaq[k]*QD0*PetscPowScalar((Vm/Vm0),ld_betaq[k]); 934 935 /* Load currents */ 936 IDr = (PD*Vr + QD*Vi)/Vm2; 937 IDi = (-QD*Vr + PD*Vi)/Vm2; 938 939 farr[offsetbus] += IDi; 940 farr[offsetbus+1] += IDr; 941 } 942 } 943 } 944 } 945 946 PetscCall(VecRestoreArrayRead(localX,&xarr)); 947 PetscCall(VecRestoreArray(localF,&farr)); 948 PetscCall(DMRestoreLocalVector(networkdm,&localX)); 949 950 PetscCall(DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F)); 951 PetscCall(DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F)); 952 PetscCall(DMRestoreLocalVector(networkdm,&localF)); 953 PetscFunctionReturn(0); 954 } 955 956 int main(int argc,char ** argv) 957 { 958 PetscErrorCode ierr; 959 PetscInt i,j,*edgelist= NULL,eStart,eEnd,vStart,vEnd; 960 PetscInt genj,excj,loadj,componentkey[5]; 961 PetscInt nc = 1; /* No. of copies (default = 1) */ 962 PetscMPIInt size,rank; 963 Vec X,F_alg; 964 TS ts; 965 SNES snes_alg,snes; 966 Bus *bus; 967 Branch *branch; 968 Gen *gen; 969 Exc *exc; 970 Load *load; 971 DM networkdm; 972 #if defined(PETSC_USE_LOG) 973 PetscLogStage stage1; 974 #endif 975 Userctx user; 976 KSP ksp; 977 PC pc; 978 PetscInt numEdges = 0; 979 980 PetscCall(PetscInitialize(&argc,&argv,"ex9busnetworkops",help)); 981 PetscCall(PetscOptionsGetInt(NULL,NULL,"-nc",&nc,NULL)); 982 PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD,&size)); 983 PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&rank)); 984 985 /* Read initial voltage vector and Ybus */ 986 if (rank == 0) { 987 PetscCall(read_data(nc,&gen,&exc,&load,&bus,&branch,&edgelist)); 988 } 989 990 PetscCall(DMNetworkCreate(PETSC_COMM_WORLD,&networkdm)); 991 PetscCall(DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(Branch),&componentkey[0])); 992 PetscCall(DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(Bus),&componentkey[1])); 993 PetscCall(DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(Gen),&componentkey[2])); 994 PetscCall(DMNetworkRegisterComponent(networkdm,"excstruct",sizeof(Exc),&componentkey[3])); 995 PetscCall(DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(Load),&componentkey[4])); 996 997 PetscCall(PetscLogStageRegister("Create network",&stage1)); 998 PetscCall(PetscLogStagePush(stage1)); 999 1000 /* Set local number of edges and edge connectivity */ 1001 if (rank == 0) numEdges = NBRANCH*nc+(nc-1); 1002 PetscCall(DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1)); 1003 PetscCall(DMNetworkAddSubnetwork(networkdm,NULL,numEdges,edgelist,NULL)); 1004 1005 /* Set up the network layout */ 1006 PetscCall(DMNetworkLayoutSetUp(networkdm)); 1007 1008 if (rank == 0) { 1009 PetscCall(PetscFree(edgelist)); 1010 } 1011 1012 /* Add network components (physical parameters of nodes and branches) and number of variables */ 1013 if (rank == 0) { 1014 PetscCall(DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd)); 1015 genj=0; loadj=0; excj=0; 1016 for (i = eStart; i < eEnd; i++) { 1017 PetscCall(DMNetworkAddComponent(networkdm,i,componentkey[0],&branch[i-eStart],0)); 1018 } 1019 1020 PetscCall(DMNetworkGetVertexRange(networkdm,&vStart,&vEnd)); 1021 1022 for (i = vStart; i < vEnd; i++) { 1023 PetscCall(DMNetworkAddComponent(networkdm,i,componentkey[1],&bus[i-vStart],2)); 1024 if (bus[i-vStart].nofgen) { 1025 for (j = 0; j < bus[i-vStart].nofgen; j++) { 1026 /* Add generator */ 1027 PetscCall(DMNetworkAddComponent(networkdm,i,componentkey[2],&gen[genj++],6)); 1028 /* Add exciter */ 1029 PetscCall(DMNetworkAddComponent(networkdm,i,componentkey[3],&exc[excj++],3)); 1030 } 1031 } 1032 if (bus[i-vStart].nofload) { 1033 for (j=0; j < bus[i-vStart].nofload; j++) { 1034 PetscCall(DMNetworkAddComponent(networkdm,i,componentkey[4],&load[loadj++],0)); 1035 } 1036 } 1037 } 1038 } 1039 1040 PetscCall(DMSetUp(networkdm)); 1041 1042 if (rank == 0) { 1043 PetscCall(PetscFree5(bus,gen,load,branch,exc)); 1044 } 1045 1046 /* for parallel options: Network partitioning and distribution of data */ 1047 if (size > 1) { 1048 PetscCall(DMNetworkDistribute(&networkdm,0)); 1049 } 1050 PetscCall(PetscLogStagePop()); 1051 1052 PetscCall(DMCreateGlobalVector(networkdm,&X)); 1053 1054 PetscCall(SetInitialGuess(networkdm,X)); 1055 1056 /* Options for fault simulation */ 1057 ierr = PetscOptionsBegin(PETSC_COMM_WORLD,NULL,"Transient stability fault options","");PetscCall(ierr); 1058 user.tfaulton = 0.02; 1059 user.tfaultoff = 0.05; 1060 user.Rfault = 0.0001; 1061 user.faultbus = 8; 1062 PetscCall(PetscOptionsReal("-tfaulton","","",user.tfaulton,&user.tfaulton,NULL)); 1063 PetscCall(PetscOptionsReal("-tfaultoff","","",user.tfaultoff,&user.tfaultoff,NULL)); 1064 PetscCall(PetscOptionsInt("-faultbus","","",user.faultbus,&user.faultbus,NULL)); 1065 user.t0 = 0.0; 1066 user.tmax = 0.1; 1067 PetscCall(PetscOptionsReal("-t0","","",user.t0,&user.t0,NULL)); 1068 PetscCall(PetscOptionsReal("-tmax","","",user.tmax,&user.tmax,NULL)); 1069 1070 PetscCall(PetscMalloc1(18*nc,&user.ybusfault)); 1071 for (i = 0; i < 18*nc; i++) { 1072 user.ybusfault[i] = 0; 1073 } 1074 user.ybusfault[user.faultbus*2+1] = 1/user.Rfault; 1075 ierr = PetscOptionsEnd();PetscCall(ierr); 1076 1077 /* Setup TS solver */ 1078 /*--------------------------------------------------------*/ 1079 PetscCall(TSCreate(PETSC_COMM_WORLD,&ts)); 1080 PetscCall(TSSetDM(ts,(DM)networkdm)); 1081 PetscCall(TSSetType(ts,TSCN)); 1082 1083 PetscCall(TSGetSNES(ts,&snes)); 1084 PetscCall(SNESGetKSP(snes,&ksp)); 1085 PetscCall(KSPGetPC(ksp,&pc)); 1086 PetscCall(PCSetType(pc,PCBJACOBI)); 1087 1088 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&user)); 1089 PetscCall(TSSetMaxTime(ts,user.tfaulton)); 1090 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 1091 PetscCall(TSSetTimeStep(ts,0.01)); 1092 PetscCall(TSSetFromOptions(ts)); 1093 1094 /*user.alg_flg = PETSC_TRUE is the period when fault exists. We add fault admittance to Ybus matrix. 1095 eg, fault bus is 8. Y88(new)=Y88(old)+Yfault. */ 1096 user.alg_flg = PETSC_FALSE; 1097 1098 /* Prefault period */ 1099 if (rank == 0) { 1100 PetscCall(PetscPrintf(PETSC_COMM_SELF,"... (1) Prefault period ... \n")); 1101 } 1102 1103 PetscCall(TSSetSolution(ts,X)); 1104 PetscCall(TSSetUp(ts)); 1105 PetscCall(TSSolve(ts,X)); 1106 1107 /* Create the nonlinear solver for solving the algebraic system */ 1108 PetscCall(VecDuplicate(X,&F_alg)); 1109 1110 PetscCall(SNESCreate(PETSC_COMM_WORLD,&snes_alg)); 1111 PetscCall(SNESSetDM(snes_alg,(DM)networkdm)); 1112 PetscCall(SNESSetFunction(snes_alg,F_alg,AlgFunction,&user)); 1113 PetscCall(SNESSetOptionsPrefix(snes_alg,"alg_")); 1114 PetscCall(SNESSetFromOptions(snes_alg)); 1115 1116 /* Apply disturbance - resistive fault at user.faultbus */ 1117 /* This is done by adding shunt conductance to the diagonal location 1118 in the Ybus matrix */ 1119 user.alg_flg = PETSC_TRUE; 1120 1121 /* Solve the algebraic equations */ 1122 if (rank == 0) { 1123 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n... (2) Apply disturbance, solve algebraic equations ... \n")); 1124 } 1125 PetscCall(SNESSolve(snes_alg,NULL,X)); 1126 1127 /* Disturbance period */ 1128 PetscCall(TSSetTime(ts,user.tfaulton)); 1129 PetscCall(TSSetMaxTime(ts,user.tfaultoff)); 1130 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 1131 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&user)); 1132 1133 user.alg_flg = PETSC_TRUE; 1134 if (rank == 0) { 1135 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n... (3) Disturbance period ... \n")); 1136 } 1137 PetscCall(TSSolve(ts,X)); 1138 1139 /* Remove the fault */ 1140 PetscCall(SNESSetFunction(snes_alg,F_alg,AlgFunction,&user)); 1141 1142 user.alg_flg = PETSC_FALSE; 1143 /* Solve the algebraic equations */ 1144 if (rank == 0) { 1145 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n... (4) Remove fault, solve algebraic equations ... \n")); 1146 } 1147 PetscCall(SNESSolve(snes_alg,NULL,X)); 1148 PetscCall(SNESDestroy(&snes_alg)); 1149 1150 /* Post-disturbance period */ 1151 PetscCall(TSSetTime(ts,user.tfaultoff)); 1152 PetscCall(TSSetMaxTime(ts,user.tmax)); 1153 PetscCall(TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER)); 1154 PetscCall(TSSetIFunction(ts,NULL,(TSIFunction) FormIFunction,&user)); 1155 1156 user.alg_flg = PETSC_FALSE; 1157 if (rank == 0) { 1158 PetscCall(PetscPrintf(PETSC_COMM_SELF,"\n... (5) Post-disturbance period ... \n")); 1159 } 1160 PetscCall(TSSolve(ts,X)); 1161 1162 PetscCall(PetscFree(user.ybusfault)); 1163 PetscCall(VecDestroy(&F_alg)); 1164 PetscCall(VecDestroy(&X)); 1165 PetscCall(DMDestroy(&networkdm)); 1166 PetscCall(TSDestroy(&ts)); 1167 PetscCall(PetscFinalize()); 1168 return 0; 1169 } 1170 1171 /*TEST 1172 1173 build: 1174 requires: double !complex !defined(PETSC_USE_64BIT_INDICES) 1175 1176 test: 1177 args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason 1178 localrunfiles: X.bin Ybus.bin ex9busnetworkops 1179 1180 test: 1181 suffix: 2 1182 nsize: 2 1183 args: -ts_monitor -snes_converged_reason -alg_snes_converged_reason 1184 localrunfiles: X.bin Ybus.bin ex9busnetworkops 1185 1186 TEST*/ 1187