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